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SUMMARY 


MSFC is participating in the current U.S. resurgence of interest in new 
liquid propellant rocket propulsion systems. In the course of this work 
modification of existing MSFC test stands is in process. As a part of the 
activation of the new liquid propulsion test facilities it is necessary to analyze 
total propulsion system stability. To accomplish this, several codes have been 
built to run on desktop 386 machines. These codes enable one to analyze the 
stability questions associated with the propellant feed systems. In addition, the 
work of Dr. Mitchell at Colorado State has been adapted to this computing 
environment and furnished along with the other codes. This latter inclusion 
furnishes those interested in high frequency oscillatory combustion behavior 
(that does not couple to the feed system) a set of code for investigations of 
proposed liquid rocket engines. 
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INTRODUCTION 


As the decade of the 1980's ended there was a resurgence of interest in 
the United States in advanced liquid propulsion systems. This resurgence 
manifested itself in the inauguration of a number of programs by both the USAF 
and the NASA. These programs take the form of both analytical work to take 
advantage of the vast increases in computer power which occurred during the 
decade of the 80's and some experimental work with new or modified engines. 
At the Marshall Space Flight Center (MSFC/NASA) test stands which have lain 
dormant since the sixties are being activated and modified. Because of the 
extensive modification of the test stands e.g., new tanks and lines, and the new 
engine design(s) to come, the MSFC needs the capability to predict stability of 
the new stand and engine combination(s). In addition, as research continued it 
became apparent that several codes for intrinsic engine stability (i.e., the high 
frequency) investigations were available and could be adapted for execution on 
the same PC 386 class of machine as the other stability codes developed under 
this effort. Thus the thrust of this work has been to gather the specifics of the 
technology applying to all of the major anticipated stability problems attributed 
directly to a propulsion system that might occur on the new or revamped MSFC 
test stands, install them on 386 class desktop computers and apply them (as 
first users) to the new stands and engines taking form at MSFC. 

References in this field are legion. Text books such as references one 
and two are available. A well known comprehensive coverage of the work 
available through the early 1970's edited by Messrs. D. T. Harrje and F. H. 
Reardon, is found in NASA SP-194 (reference three which contain 778 
references). Discussions with engine manufacturers and with Dr. Reardon of 
the University of California in Sacramento reveal a consensus that stability 
prediction theory has not changed a great deal since SP-194 was published. 
Thus the low and intermediate frequency stability analyses are based upon that 
technology. We have drawn upon the very recent work of Professor Mitchell at 
Colorado State (USAF funded; see references four, five and six) for the high 
frequency stability modeling and prediction codes. 
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In his two volume report, Professor Mitchell and his group have extended 
greatly their earlier work (see reference four). This report contains 54 
references, a fair number of which are dated after those of SP-194 i.e., after 
circa 1968. In this latest work, combustion chambers of rectangular as well as 
circular cross section, are analyzed. In addition, quarter wave, Helmholtz and 
arbitrary geometry (L,T etc.) dampers are covered as are damping chamber 
liners. In this latest stability code virtually any combination of these devices may 
be present simultaneously in the combustion chamber in rather arbitrary 
locations e.g., there may be a partial liner in the combustion chamber while 
there are simultaneously resonators around the periphery of the injector end of 
the combustion chamber. Because Mitchell's work covers only high frequency 
stability questions dealing with combustion chamber design and not test stand 
design or modification it is of most interest to engine designers or those 
monitoring the development of new engine designs (or the modification of 
existing engine designs). 

Professor Mitchell and the USAF (Mr. J. Levine AFSC Edwards Air Force 
Base) have been most helpful in providing access to his work especially by 
providing copies of his source code on computer media. 

Our primary goal then was to create a family of programs running on 
tabletop computers that are of immediate use to engineers designing test 
stands destined for use with new liquid propellant propulsion systems. These 
programs have been formulated and installed in such a way that engineers 
without exhaustive backgrounds in fluid dynamics, stability theory and 
advanced mathematics will find them useful to the point of turning to and using 
them during the design and the checkout phases of their new propulsion 
hardware initiatives. In addition, we have concatenated to these programs 
those of Professor Mitchell, modified for the desktop environment, for the use of 
those people interested in high frequency combustion stability questions. 
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THEORY 


Directly associated with an engine and propellant feed system there are 
generally taken to be three separate stability regimes. They are classified by 
the frequencies at which they tend to oscillate. Reasonably enough they are 
called the low, intermediate and high frequency regimes. 

While in concept the analysis of all three regimes could be lumped into 
one extremely large set of code, in practice separate modules or programs are 
used for each regime. A further distinction between the regimes occurs in that, 
as will be developed, the low and intermediate regimes involve an intrinsic 
interaction between the combustion chamber processes and the propellant 
feedlines while the high frequency regime phenomena are usually confined to 
the immediate environs of the combustion chamber i.e. the combuster. 

For these reasons the test stand designers may well want to pay 
particular attention to the two lower frequency stability regimes. 

In what follows linear analysis is used. This is permissible if the 
quantities in question (e.g. mass flow rate and pressure) are assumed not to 
deviate very far from their set point (mean or average) values. Thus stability is 
evaluated, for example by investigating whether or not a pressure or mass flow 
rate, once perturbed, will return to its average or mean value i.e., whether or not 
the perturbation or small deviation vanishes with time. Because of this linear 
approach all the analytical tools of linear stability analyses are available. Also, 
because of this approach a number of analyses of any given hardware system 
will have to be performed; each corresponding to a different operating condition 
of the propulsion system. Different conditions would include such things as 
different amounts of liquid in the propellant tanks, operation with different 
mixture ratios and different valve settings. 
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THEORY OF PROPELLANT FEED LINES 


To characterize dynamically a propellant feed system it is necessary to 
treat it in the same vein as other circuit problems. The circuit elements 
encountered may be categorized as follows. Distributed lines (incorporating the 
effects of distributed hydraulic capacitance and inertance), discrete hydraulic 
capacitors (used to model the engine manifolds and the propellant tanks), 
hydraulic resistors (which are used to characterize such things as the injectors 
through which propellant is introduced into the combustion chamber) and 
propellant pumps. 

The methodology suggested in section 5.4.1 of SP-194 is used as 
follows. The admittance, G, at any point in the feed system is defined to be, in 
general, a complex quantity given as the ratio of the mass flow rate perturbation 
to the local pressure perturbation. Also, in general, any "downstream" 
admittance is a function of an "upstream" admittance. To illustrate the last point, 
it will be shown that the admittance of the feedline attached to a tank is a 
function of the admittance "looking into" the tank. With the foregoing in mind, 
the admittance as seen "looking" from the combustion chamber system toward 
the tank is formed by the following product of admittance ratios. 



where each ratio represents the admittance ratio associated with a given 
feedline circuit element. Typically one starts at the tank, considered to be a 
large, here constant pressure, liquid source and work one's way down the 
piping to and through the engine injectors. These admittance expressions are a 
key element in evaluating both the low and intermediate stability modes of the 
engine - propellant feed system. 

In the work which follows immediately below, it is convenient to have at 
hand a general expression for the input admittance ratio of a certain network. 
This expression may be used to characterize all the combinations of discrete 
hydraulic circuit elements of interest. This general network is shown below 
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All the circuit elements will be taken to have sign convention 

corresponding to passive circuit elements. As will be pointed out below this 

JL 

needs very special interpretation in the case of a pump. Combining Z2 and Gi 
in parallel and then summing with A (after a little algebra) produces 



Note that if it is desired to open up a circuit element (i.e., no flow through it) then 
the appropriate Z is set equal to 00 (or G to 0 ). If a circuit element is to present 
no obstacle to flow then its Z is set equal to 0 (or G to°°). 

HYDRAULIC RESISTANCE 

The first derivation will concern itself with an energy dissipating hydraulic 
resistance used to characterize an orifice. First consider the development of the 
linearized relationship between pressure drop and liquid flow through an 
orifice. This relationship is usually characterized by a square law relationship, 
i.e. 

p= Rm 2 

such that 
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The incremental resistance to flow at the operating point (defined by P 
and m )is taken to be the slope of the p versus rh curve evaluated at P=P and 
rh=m. 

Thus 


r *= 


dp ip=p 
drh m=m 


and from before 



Substituting yields 


r’*=^-(Rrii 2 ) = 2Rrtil R(p 4- m) 
aril m=m 


so 


r'= 



In all this work the pressure is normalized by the average combustion pressure 
(Pc) and the mass flow by the total p ropellant average mass flow Thus the 
normalized resistance is 
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r’= 2=2-1 [m 


mj \PC 


where it<L is the average mass flow rate through the resistor. This derivation 
holds for a short orifice i.e., one in which the transit time of the fluid through the 
orifice is short. With this expression for resistance in hand, the admittance ratio 
is developed. A diagram is helpful at this point. 


Zo 

o AMA o 


o 


Orifice Circuit 



From the figure it is seen that z i-Zo and Z2-°° j n the general formulation. Thus 
the input admittance ratio is 

G2- 1 

Gi 1+ ZoGi 


HYDRAULIC CAPACITANCE 

Hydraulic capacitance is needed to account for the compressibility of the 
propellant liquids. For a lumped capacitance an analogous development to that 
of the resistor holds with the exception that time rates of change of pressure 
must be introduced i.e., the dynamics of fluid compression govern the dynamic 
behavior of the element. By definition the dynamics of a lumped capacitor are 
governed by 

dp'_ m’ 
dt c 

Under the assumption that forced solutions are of paramount interest here, the 
initial conditions applying to P' and m' are assumed zero yielding the LaPlace 
transformed equation as 


P(s)= 


Ift'(s) 
c s 


8 




where 


s = LaPlace Operator 
p’(s)=L[p’(t)] assumed to exist 
m’(s)=L[m'(t)] assumed to exist 
Once again a picture is helpful as 


o 


o 



Hydraulic Capacitor Circuit 


As before, from inspection, 


Z 1= 0 and Z 2 =^ 


substitution yields 


Cj2_ 

Gi 


1+£L 
Gi . 


Two observations are in order. First consider the defining differential equation 

dp'_ m' 
dt c . 


From this it is evident that for c very large, as might occur with a propellant tank 
characterization, it will be difficult for p' to change. A second observation is 
evident from the transform of the differential equation i.e., 


p’(s)= 


ifa'(s) 
cs . 


Letting s=jvv to consider frequency response effects shows 


p'(jO))= 


m'(jco) 

jcoc 
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Thus as Q)-> 0 when characterizing a tank then care must be taken if 

the frequency of interest is very low. In that case additional characterization 
should be considered e.g. perhaps an orifice or other purely real admittance 
should be added to the description. 

PROPELLANT FEED PUMP 

Following the methodology of SP-194 a lumped parameter 
characterization of a pump may be formulated. Some contemplation of a pump 
is first in order. The same pump may, depending upon the operating point, 
have either a rising, flat, or falling pressure change versus fluid flow 
characteristic. Graphically this may be shown as below. 



From the figure it is clear that the slope of 

dp' 

dih' 


may be greater, equal to or less than zero (this sign may very well vary with a 
given pump's operating point). Thus the pump impedance corresponding to the 
head-capacity relationship may be taken to be a resistor but one whose sign 
may not be "taken for granted". Clearly there is fluid in the pump. Therefore, it 
is necessary to associate with a dynamic pump model some inertance and 
capacitance. Numerical values of these quantities probably have to be gotten 
from dynamic testing. Such things as cavitation effects could cause values all 
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out ol relationship to the nominal amount of fluid in the pump. The equivalent 
pump circuit is as shown. 


Zp L 

o — AAAA-® 


q 



Gi 


Hydraulic Pump Circuit 


By inspection 


Zi-Ls+Zpi Z2— 


thus 


G2_ . 


1 ■ cs 
Gi 


1+F- 

Gi 


Gi 1 +{Zp+Ls)cs+(Zp+Ls)G 1 l-KLs+Zp)(cs+G,) 


or 


1 1 cs 

G2_ Gj . 

Oi l+Z p G 1 |l+^-)(l+tps)' p Zp 


Note: in this formulation (to be consistent with the general formulation) 



It may very well be that 


Zp= 


dp |P=P_ 
9 rfa 


0 
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In which case a negative number should be inserted for in the admittance 
ratio. 

PROPELLANT DUCT 

Another major non-active propellant delivery circuit element for which the 
admittance ratio is needed is the pipe or duct which carries the liquid from 
element to element. The following development is for a pipe of uniform cross 
section. Pipes of nonuniform cross section may well require the application of 
an ad hoc numerical solution for each one (there are some exceptions e.g., 
incompressible flow in a rigid wall pipe). The following development will be 
started including in the development the effect of a resistive pressure drop. 
However initial calculations on one proposed MSFC test stand design showed 
very low values of such drop and hence at the appropriate time they were 
dropped as being negligible. However the development can be readily 
expanded to include the resistive effect. 

Partial differential equations may be used (for an alternate treatment see 
Appendix A) to treat the distributed nature of the feedline. One could, of course, 
approximate the feedlines by an interconnected set of discrete resistors, 
capacitors and inertances. But aside from the approximation (which could be 
made as negligible as desired by adding more elements) it is arguably more 
computationally efficient to use the results of the continuous representation. 

The dynamic characterization of a constant cross section fluid duct with 
rigid walls will be accomplished with three objectives in mind. The first will be to 
characterize a lossless duct to obtain its admittance ratio. Then a short 
exposition will be included to indicate the straightforward way in which the 
analysis may be expanded to include dissipating (distributed resistance) effects. 
Lastly will be developed expressions for the standing pressure and flow waves 
to be found along a duct (pipe). 

The dynamics of the pressure flow relations in a duct are governed by the 
one dimensional wave equation. The development is begun by deriving the 
wave equation. Refer to the figure which depicts an infinitesimally small 
segment of a uniform duct 
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where 

A 

Z=per unit length series impedance (resistance, inertance) 

A 

y=per unit length shunt admittance (typically conductance, 
capacitors) 

note that x is measured from the receiving end of the duct. 

By inspection the pressure drop across the element, along the duct is 

dP=mzdx or^=mz 
dx 

and as the liquid flows through the element there is a decrease 

drh=Pydx or^=yP 
dx 

Note in these equations that in general P -^-y- anc * z are either functions of x or 
the LaPlace variable S or both. The most usually used expression for Z in 
transmission line network theory is probably 

z=Ls+R 

and for the admittance 

y=G+cs 
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Some discussion of y is in order. The use of shunt capacity and 
conductance is in accord with usual network theory models. However, unless 
the duct leaked uniformly along its length, G, the conductance, would be set 
equal to zero for the purposes of this work. But it should be borne in mind that 
the use of other admittance functions has been suggested. For example, if a 
duct has internal damping such as might be provided by a liner which could be 
modeled as a set of distributed Helmholtz resonators then the following model 
has been proposal (in SP-194 on p. 274, Section 6. 2. 5. 4) for the shunt 
admittance. 



In this case y may be computed as follows. 



G 2 => y 


y_ LnC n S 2 +R n CnS+1 | Q bS 
C n s 
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(L n C n +CbC n )s^+RnC n s+1 

y= CnS 

Returning to the overall derivation the derivative of each equation is taken with 
respect to x and the results substituted to obtain the one dimensional wave 
equations for pressure and flow as 


d 2 P 

dx 2 


=zyP 


d d rh 

dx 2 


=zyrh 


Assume a solution for the pressure of the form 


P=Ae ax 


thus 


2 

= Aae ax and = Aa 2 e ax 
dx dx 2 


Substituting into the wave equation yields 
2 

= Aa 2 e ax = zyAe ax 
dx 2 

from which it is seen that 
a 2 = zy or a = ±Vzy 

Thus the traveling wave expression for pressure becomes 
p = Ae^r x + B e-^ x . 


Recalling the original expression 


15 



m = 


idE 
Z dx 


and substituting for P the expression just obtained yields 
m = e^ x - e-^ x 


The constants A and B are evaluated from the end conditions. Let p r and 
be the pressure and flow rate at the receiving end of the duct. Then recalling 
that here x is measured from the receiving and let x = 0 to obtain 


P r = A + B 

m r = Ayj -B^J 


solve simultaneously for A and B 



Substituting into the expression for pressure gives 
p = + ,j, r yz[ 

P = Pr cosh xvz y + sinh xvzy 

Similarly for the mass flow rate 

m = m r cosh xVzy + sinh xVzy 

These equations have been written when the distance is measured from 
the receiving end. These equations can be written in terms of sending end 
pressure and mass flow rate and with distance measured from the sending end 
provided -x is substituted for x in the equations above. 
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Because in the above equations positive x is measured back along the 
duct from the receiving end (with x = 0 at the receiving end), distances along the 
line (with x = 0 at the sending end), are to be regarded as -x. Recalling that 
cosh is an even function and sinh is an odd function then substituting -x for x 
yields 

P = P s cosh xVzy - sinh xVzy 

m = ms cosh x-fzy - P s / )j^ sinh xVzy 

for a sending point |x| distance from the sending end. It is convenient to define 
the surge or characteristic line impedance as 

Zo= #. 

Letting the length of the duct be 1 and the ratio of Pr to ihr be Zr then the input 
impedance of the duct as seen from the sending end is the quotient of P s an d * s . 

P s = P r cosh {Vzy +rhrZo sinhlVzy 

m s = rh r coshlVzy +z-i sinh IVzy 
Zo 

z _Ps__ z Z r cosh IVzy +Z 0 sinh IVzy 
s ms ° Zo cosh lVzy +Z r sinh IVzy. 

To obtain the result in SP-194 (5.4.1 -2) divide by coshlVzy 

Z r + Zp tanh jVzy 
Zo + Z r tanh IVzy 

and factor out Zo and divide by Z r 

ZlL 

4 


1 + ^P-tanh IVzy 
Zr 

1 +^tanh IVzy 
Zo 
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Let 

Gi = 4" and G 2 = j- 
Z r Z 2 

then 

Zs _ Gi 
4 G 2 
or 

1 +^-tanh l-/zy 1 + — 1 tanh l-/zy 
G 2 _ Zp G-|Zp 

Q 1 1+^0-tanhlVzy 1 + GiZ 0 tanh IVzy 
Zr 

This is the desired admittance ratio which allows the duct to be treated in the 
same manner as a discrete element in the propellant teed line. 

In modeling one proposed MSFC Test Stand piping layout a number of elbows 
were encountered. It has been shown by Jackson (as reported in SP-194, 
Section 3.2.2) that the inertance of an elbow can be regarded as the sum of two 
contributions. One of these is the inertance associated with the fluid treated as 
if it were flowing in a straight duct. The other is an increase in inertance caused 
by the duct’s curvature. The analytical and experimental results of Jackson are 
shown as a curve of the latter contribution as a function of duct diameter and 
parameterized by the angle through which the duct is bent. The curve is 
reproduced here for convenience. 
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The approach taken by this author is to derive an equivalent lumped series 
circuit element. This element models the curved duct as a straight pipe with a 
decreased area and an increased length compared to a straight pipe having the 
length of the curved duct's center line and having a decreased cross section 
area compared to the bent duct. In this development a constant volume of the 
two ducts' fluid is maintained. The deviation is a follows. Note that the 
expressions for hydraulic inertance and capacitance for a lumped line are given 
by 

L= 1/A 

c = (2) 1A 


Thus for two straight ducts of the same volume one may write 
expressions for inertance and capacitance for each as 


C, = ( 2 ) 1,A, C 2 = ( 2 ) 1 2 a 2 


Ai 



For a constant (lumped) volume of fluid the compressibility will be 
the same regardless of its form. Therefore, 


Ci = C2 


and thus 


liflj = l 2 fla 


or 


k =0'- ■ 
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Now let 


L2 = XLi 


Where X is defined by this equation. Then substituting m in the inertance 
relationship yields 


h 

A2 



from which 

^fr 2 


from before 
substituting 



thus the "new area” of the equivalent line becomes 



and the "new" line length of the equivalent line is 


All that remains is to evaluate 



X 


. This is done by recalling 


where 


L '=f 

Ai 


; Lj + L = L2 


and L" is found from the graph. Our current experience indicates that because 
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the effect of the added inertance is small especially at low frequencies it speeds 
computation to ignore this effect upon first analysis and then add its effect at the 
high frequencies at a later date to evaluate its effect upon stability. 

FEEDLINE PRESSURE AND FLOW RELATIONSHIPS 


During the course of this work the question of where to place pressure taps in 
the feedlines arose. This question came about because previous experience 
showed the importance of location when measuring pressure in a feedline. 
Fortunately all the methodology developed to this point is applicable. Recall 
that expressions for the admittance (and hence impedance) at any circuit 
element location (hydraulic capacitor, inertance, valve of duct end) "looking 
back" toward the tank have been developed. Similar methodology could be 
used to develop the admittance of impedance expressions "looking 
downstream" toward the engine's combuster. With these recollections in mind it 
is seen that it is possible to fashion a model of a duct being driven by an 
excitation source in series with a source or generation impedance and 
terminated in a load or receiving end impedance. 

In the interest of additional pedagogy and to develop additional notation the 
pressure flow relations applicable to a duct will be rederived. However, the 
derivation will be done in such a way as to expose the functional relationships 
between pressure and flow and such things as excitation description, source 
and load impedances and the line parameters. 


From before the duct is treated as a transmission line. 



RAx 


AW 


GAx 


a 


CAx 


LAx 

W 


dm 

Q 


P + 


dp 

Bx 


Ax 



21 



Recall the previous discussion about the various impedances e.g. in cases of 
smooth ducts it is expected that G ->o, R is small and the dynamics are 
dominated by L and C). If other conditions obtain then the methodology 
followed here may be used to solve the other case. Assuming the first term of 
the Taylor's Series expansion is adequate to describe the change in pressure 
and mass flow rate the differential equations for this infinitesimal length of duct 
are formulated by equating the source pressure (at x) to the pressure drop 
across the duct element summed with the pressure at the receiving end of the 
duct. 

p(x,t)= RAxm(x,t) +LAx -°^ X,t ^ + p(x,t) + — fo-’^ Ax 

ax. ox. 

Likewise the equation for the flow rate relationship is formulated by 
equating the input flow rate to the element to that flowing through the shunt 
admittance summed with the flow leaving the element. 

m(x,t)=p(x,t)GAx + cAx ^PQ ^ ) + m(x,t)+ A x - ^P^ x,t ^ 

dt dx 


From these two equations it follows that 

^a=Rri 1( x,. )+ ii5#a 

dx dt 


-dm(x,t) _ 


dx 


= Gp(x,t) + c 


dp(x,t) 

dt 


are the defining equations. To transform from the time domain to the frequency 
domain (and hence transform from differential equations in time to algebraic 
equations in the complex frequency variable s) the La Place transform is 
utilized. Initial conditions will be set equal to zero on heuristic grounds. Thus 
the effects of driven or particular solutions may be studied as the problem is 
formulated. If a given set of initial conditions were known the equations may be 
rederived to include them. The LaPlace transformed equations are 

dP(x,s) _ ^ + Lsm(x,s) = (Ls + R)m(x,s) 

dx 

-dm( x , _ s) = Gp(x>s) + csP(XjS ) = (cs + Q)P(x,s) 


Taking the partial derivative with respect to x of each of these equations and 
substituting the expression above for the resulting first partials yields. 

a2p(x,s) - (R + Ls)(G+cs)P(x,s)=0 
dx 2 
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d 2m( * ,s) - (R + Ls)(G+cs)m(x,s)=0 
9x 2 


now define 

n 2 = (Ls + R)(G + cs) : the propagation constant and note that for a lossless line 
(a usual first approximation) 
n 2 = Lcs 2 


To achieve a solution, one is assumed, analogously with previous work. 
m(x,s) = Ai(s)e _nx + Bi(s)e nx 
P(x,s) = A 2 (s)e- nx + B 2 (s)e nx 

Note the functional arguments of the "constants." The "constants" (A's and B's) 
are evaluated from boundary conditions. Therefore, to evaluate them some 
specific configuration must be specified. 

The configuration is chosen as shown below. 


m(d,s) 



By inspection 

P(0,s) = P g (s) - Z g rii(0,s) = A 2 (s) + B 2 (s) 

P(d,s) = Zt(s) m(d,s) 

Differentiate the assumed pressure solution with respect to x 

= -nA 2 (s)e* nx + nB 2 (s)e nx 
c/x 
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and substitute it into the original transformed PDE relating the rate of change of 
pressure to the mass flow rate to yield 

m(x,s) = — l - [ nA 2 (s)e- n * - nB 2 (s)e nx ] 

(Ls + R) 


Comparing this expression with the assumed solution for the mass flow rate 
produces the relationships among the "constants" 


Ai(s) = 


nA2(s) 
(Ls + R) 


Bi(s) = 


-nB 2 (s) 
(Ls + R) 


These expressions for A i(s) anc j Bi(s) ma y p e substituted into the originally 
assumed solutions to yield 


rirfx sl = _2M!Le.nx . ”^00 
(Ls + R) (Ls + R)' 

P(x,s) = A 2 (s)e nx + B 2 (s)e nx 


*nx 


Note that 


n 


(Ls + R) 


(Ls + R)(Cs + G) ^y (Cs + GJ~ _ l 


(Ls + R) 2 V (Ls + R) A, 


where Zo is defined to be the surge or characteristic impedance and 




for a lossless line. 


using the definitions above produces 

m(x,s) = - ^|^e nx 

Zo Zq 

P(x,s) = A 2 (s)e- nx + B 2 (s)e nx 

The sending end boundary conditions are applied to yield 

P(0,s) = Pg(s) - Z g (s)m(0,s) = Pg(s) - |^A 2 (s) - B 2 (s)] = A 2 (s) + B 2 (s) 

A) 

and the receiving end boundary conditions yield 
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P(s,d) = Ztm(s,d) = |KA 2 (s)e-nd - B 2 (s)e" d ] 
A> 


These two simultaneous linear equations in the unknowns A 2 (s) 
B 2 (s) may be solved to yield 

z^- p s (s) 

A l( s)= 

1 + MNe- 2nd 


and 


( -rA-)Ne- 2nd P g (s) 


B 2 (s)=- 


Zo +z. 


1 - MNe- 2nd 


where the reflection operators are defined by 


N= Z ° 


-z. 


Zq + Zt 


M = 


Zq - Zfi 


Zq + Zg 

Using these solutions for A 2 and B 2 and substituting them into the assumed 
solutions for pressure and flow rate produces the transfer functions between the 
pressure generator p g( s ) and the pressure and flow perturbations at a given 
distance down the duct as 


P(x,s) Zp+Zj 

P«(s) ' 


_[e-nx . (|olA )e -n(2d-x) ] 


Z 0 +Z/ 


[1 . f (Zq - Zt) (Zq - Zg) , 

1 l (Zo+Zt)(Z 0 +Z g ) J J 


and 


M(x,s) _ Z o +Zj 
P g (s) “ 


■J— [e-^ + (|^A)e^(2 d ^)] 


'Z 0 +Z t / 


_ (Z 0 -Z t )(Z 0 -Z g ) 

1 l (Z o +Zt)(Z 0 +Z g ) ;C 1 


While there are special cases of interest (some considered below) it is seen, 
especially recalling that all the z's are functions of s (s = oc + jco ) t that the 
solutions following from this transfer function, given a generating function are 
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generally of a complicated nature. This fact leads to numerical investigations to 
explore given cases. But before investigating numerical approaches consider 
some special cases. 

The most perfect (and impossible) case occurs when the duct is considered 
lossless and the reflection operators are zero (the matched load, lossless case 
i.e., z e = z t=Zo ). Under these assumptions the transfer functions yields 

P(x,s) = P g (s)e- nx 

This indicates that an input pressure signal will travel down the duct undistorted 
in wave shape and only delayed in time (transport lag phenomenon) or, from 
the frequency domain viewpoint, there is introduced a phase lag proportional to 
the distance down the duct. 


Another special case of some interest occurs when G = 0 and R « Ls. The 
problem to be attacked is that the quantity n contains s in a square root form. 


n = Y(Ls + R)(G + Cs) 


Under this set of assumptions 

n = Y(Ls + R) Cs = 1 4 — VLC = VLCs(l +^-) 

or 

nsVECs + lRfyC) 


Once again, assuming a matched load i.e. duct terminated in the surge or 
characteristic impedance Zq and z g = Zq ^ p ressure as a function of distance 
and the generation pressure is given by 


P(x,s) = P g (s)e-^C s e' ( 2 )R VY x 


The inverse transform to get to the time domain is obtained by inspection as 
P(t,x) =e' ( ^ )R V^ x Pg(t -VLCx)u(t -YLCx) 
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where u(t-VECx) j S the familiar unit step function notation used as a sectioning 

function so that no time response occurs before t=VLtTx i.e. this is a causal 

system (therefore physically realizable). Thus this wave may be interpreted as 

1 

an attenuated wave traveling at a speed = VLC , having zero value until 
t> VlC'x Notice that the form of this solution is the same as that 

encountered previously for the lossless line with the exception that here there is 
attenuation i.e. the wave maintains its shape as it moves along but it shrinks in 
amplitude. 

There are other combinations of the duct parameters which yield a (possibly) 
attenuated but undistorted traveling wave. These other cases are not of 
particular interest here because they do not assume G to be zero. In any case a 
duct having parameters so designed is called a "distortionless" or "heaviside" 
duct. 

Before moving away from duct characterization the subject of transient events 
will be addressed. There are a number of catastrophic events that could be 
postulated as transient events. However, the attention here will concern itself 
only with a step change to a different load. This might be used to model fairly 
rapid changes in valve settings (openings or closings). The general approach 
to handling any transient event of this generic type will always involve solving 
the given forced-solution for the first configurations' solution which in turn 
become initial conditions for the solution of the second configuration etc. To this 
point the solutions of the wave equation where these have involved transforms 
have assumed the initial condition to be zero (to get transfer functions, for 
instance). To relax this assumption return to the two first order PDE’S which 
define the pressure and mass flow rate dynamics when they are LaPlace 
transformed but this time include the initial conditions. 

= Rm(x,s) + Lsm(x,s) - Lm(x,0) 
dx 

~ 9 ” (X,S) = GP(x,s) + CsP(x,s) - CP(x,0) 
dx 

It is of course recognized that the last term in each equation is caused by the 
indicated initial condition. These equations are now manipulated as before to 
yield the wave equations in mass flow rate and pressure. 
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a2p(x ’ s) - n 2 P(x,s) = L^ X,0) - C(Ls + R)P(x,0) 
dx 2 dx 

- n 2 m(x,s) = - L(Cs + G)m(x,0) 

dx 2 ox 

Written in this form all the terms on the right hand sides (rhs) are due to the 
initial conditions. One now has a choice of how to proceed. Which way is 
chosen depends upon the particular analyst and what computational or other 
insightful factors are involved. The conceptually most straightforward way is to 
terminate one solution at the switch time, determine the pressure and mass flow 
rates in the duct , use these as the initial conditions for the new configuration 
(having written its equations) and solve another subsequent problem etc. There 
are other approaches. One of these is by the use of cancellation sources which 
model the pressure and/or flow rates of the switches after the configuration 
change. This is done by introducing ideal pressure (hence no generation 
source impedance) and ideal mass flow rate sources (hence infinite source 
impedance). This technique is especially useful for solving switching modeled 
configurations changes. It converts from a problem with initial conditions to two 
or more problems without initial conditions but whose total response is obtained 
by summing the response of all the problems. 


Returning to the original question of where to place pressure taps in the line in 
order to measure standing waves in the duct some inspection of the p(x,s)/Pg(s) 
transfer function is in order. Standing waves are associated with steady state 
sinusoidal duct responses. This type of response may be examined by 
restricting s to the imaginary axis in the s-plane i.e. let s=jw. One notes in 
passing that transient information is not lost by so doing. Making an appeal to 
the Weierstrass Approximation Theorem allows one to think of the transfer 
functions as approximated by a ratio of polynomials in s, for a given x. Provided 
the resulting transfer function is analytic in the right hand plane and on the 
imaginary axis an exact inverse Fourier Transform may be written as 


h(t) = 



(Real H(jco)cos tcojdco 


t>0 


where h (t) is the impulse response of the system and 
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L [ h (t) ] =H (s) 


Having let s=jw it is seen that the numerator of the duct transfer function consists 
of two terms of phasors which may be visualized as a magnitude and phase 
angle in the complex plane. The one term is 

e-iix 

which for the example of a lossless line would be represented by a unit phasor 
at an angle determined by the frequency of interest, the duct parameters and 
the location chosen on the duct. 

The second term is more involved but it also may be interpreted as a phasor in 
the complex plane for a given x. It is not unreasonable to expect that the 
difference of these two phasor quantities would experience maxima and 
minima. As a matter of fact for a lossless duct it may be shown that for a variety 
of conditions the pressure and mass flow are displaced along the duct by one 
quarter wavelength and both are zero at points separated by one half a 
wavelength. 

As a practical matter a specific configuration, one supplied by MSFC (a 
proposed LOX feedline), was used as an example to explore the transfer 
functions from a numerical standpoint. Two problems arise as far as visualizing 
the solution is concerned. One is that of presenting the results as a function of 
both position down the duct and as a function of frequency. The solution 
suggests itself as a three dimensional plot with transfer function magnitude, 
frequency and duct position as the three coordinate axes. Fortunately 
MATHEMATICA has a built in 3D plotting capability which has been enhanced 
in the latest release. In addition its contour plot capability proved to be 
extremely useful in determining the granularity of the solution. The second 
problem is that of the granularity. The duct chosen is 62 feet in length and 
somewhat arbitrarily a frequency range of zero to six thousand radians per 
second was chosen initially . At this point the age old question of granularity 
across the 2D grid of displacement versus frequency arose. 
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An initial range was chosen as shown with the results shown in the figure 
below. 



MAGNITUDE EATIO 



nsquncT 


LOX Feedline Pressure Transfer Function 3D Depiction 

Fig. 3D- 1 


The number of points were then increased and Fig 3D-2 produced. 



magnitude iatio 
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DISTANCE 


LOX Feedline Pressure Transfer Function 3D Depiction 

Fig. 3D-2 
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Clearly a great deal more detail has arisen. The cost is clearly in sharply 
increased computing time and, eventually, saturated memory. At this point a 
contour plot of the data was run producing figure 3D-3. 



Fig. 3D-3 
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Note that many of the "ridge lines' of response peaks are smeared together e.g. 
the ones at 3250 and 1575 radians per second while it was unclear what was 
happening below about 750 radians per second. Consider that usually there is 
more energy in lower frequency modes than higher frequency ones and also 
that one usually expects a half wave to develop between the two ends of the 
duct. At this point the frequency range was reduced to zero to one thousand 
and the same number of points plotted as before. This yielded fig. 3D-4. This 
figure was created using an updated version of MATHEMATICA which allowed 
labeling the axes. Also the plotting code is included to illustrate the use of this 
program for this purpose. 

Fig. 3D-4 

ln[l]:m 

Plot3D[Abs [ ( (3. 28+I*u*32 . 37*10*- 4) / (4 . 33+I*u*32 . 37*10 A -4) ) 

* 

(Ixp[-l*I*4 . 3*10 A -4*u*x] - ( (-1+1*2 . 76*u) / (1+I*u*2 . 76) ) *Exp[ 
- 1 * 

X*u*4 .3*10 A -4* ( 124-x) J ) / (1- ( ( (-1+X *u*2 . 76) / (1+I*u*2 . 76) ) * ( 
(2.23+I*u* .0034) / 

(4 . 33+1 *u* . 0034) ) * Exp ( - 5 . 33*10 A -2*X*u]))],(u, 0,1000), (x,0, 
62 ) , Plo-tPoints 

->30, PlotLabel->"LOX Feedline Pressure Transfer Function 
30 Depiction" 


LOX Feedline Pressure Transfer Function 3D Depiction 

Fig. 3D-4 
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In this figure the half wave of pressure becomes discernable below 150 
radians/sec and much better resolution of the other lower modes occurred. 
Therefore the upper frequency value was once again reduced, this time to 400 
and fig. 30-5 resulted. 



LOX Feedline Pressure Transfer Function 3D Depiction 

Fig. 3D- 5 


This figure clearly depicts the three lowest frequency modes in reasonable 
detail. These types of figures could be developed for any other band of 
frequencies of interest. To. develop additional information including numerical 
data contour plots were made for the frequency ranges 0-150; 200-300; and 
300-400/ In addition a series of 2D plots of amplitude versus length for the fixed 
frequencies 120, 240 and 358.3 radian per second were made. These latter 
plots were made to depict the standing pressure waves corresponding to the 
first three lowest modal or resonant frequencies. They have been grouped 
together below for ease of association and cross interpretation. Once again, for 
each type of plot one representative set of code is included. 
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ln(S]:» 

CoatourFlot [Aba [ ( (3 . 28*X«u»32 . 37*X0 A -4) / (4 . 33*I«u*32 . 37*10 
<*ip!-l*X*4 . 3*10 A -4*u*»] - t t-i + X-a . 7 4*u> / (X«Z*u*2.7C) ) ***pl 


X-J*4 .3*10*-4* <124-«> ] )/ <1- < ( (-1*1 *U* 2. 74)/ (1+I*u*2.76> ) ■ ( 
(2 . 2 3 , *X*u* .0034) / 

(4.33*X*u* .0034) ) *E«pl”5.33*10 A -2*I*u] ) ) ] » {u,0, 130) , ( « , 0 , 6 
2 ) , FlotFointa 

->30, FlotLei>el->"l.OX Feedline Freeeure Transfer Funotion 
3D Depiction" 



L0X Feedline Pressure Transfer Function Topological Depiction 

Fundamental Frequency 
Fig. 3D-6 


u«X20 

FXot[W5*(((3.2»^I'>U»32.37«X0*->4)/(4.33^X*u*32.37»X0‘‘-4))« 
(I*p(-l*I*4 . 3*10*-4*u*«] * ( ( - X ♦ X • 2 . 7 4 • u ) / (l+X a U*2. 74) ) * Kxp [ 
-X* 

X *U • 4 . 3*10 A -4* ( 1 2 4 -a ) ) ) / (X- ( ( (“XoX*u*2 . 74) / ( 1 ♦ I *u *2 .76) ) • ( 
(2 . 23*X*u* . 0034 ) / 

(4.33tX*u* .0034) ) **ap [-S.33*10 A -2*X*u]))],(x,0,62),FlotFoi 
at • 

->30, FlotL«~bel->"LOX FMdXiat Irtuur* Transfer Function 
2D Depletion*] 



L0X Feedline Pressure Transfer Function 2D Depiction 

Fundamental Frequency 
Fig. 3D-7 
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LOX FeedJine Pressure Transfer Function Topological Depiction 

First Harmonic 
Fig. 3D-8 



LOX Feedline Pressure Transfer Function 2D Depiction 

First Harmonic 
Fig. 3D-9 
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LOX Feedline Pressure Transfer Function Topological Depiction 

Second Harmonic 
Fig. 3D- 10 
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From these plots it is readily apparent that several instrumentation issues arise. 
Amongst them is the fact that instrumentation placed to measure a given mode 
may well miss another mode. For example, if it were desired to measure the 
lowest frequency modal pressure and the instrument tap is put at mid duct the 
pressure due to the next mode which nulls at mid duct would not be measured. 
Thus if it were desired to measure the next modes’ pressure taps at 15 or 45 
feet (approximately one quarter of the way down the duct) would be 
appropriate. And so forth for the other modes. Using the contour and 2D plots 
numerical data can be developed to aid in placing the taps. It also is helpful in 
specifying the bandwidth of the instrumentation system. It would seem to be 
good engineering practice to record the output of the instrumentation and 
perform a spectral analysis on it to confirm the modal frequency and to make 
sure that the desired modal pressures were measured and evaluated. 

INTRODUCTION TO COMBUSTION STABILITY 

Combustion driven or coupled stability questions associated with liquid fueled 
rocket engine systems have been amply shown to be accurately categorized as 
having a high frequency regime and a lower frequency regime. In addition the 
lower frequency regime is subdivided in two regimes, sometimes called the low 
and intermediate modes of oscillation. 

The high frequency mode of oscillation is exclusively associated with a specific 
combuster design. It depends greatly on such things as injector design, 
combuster shape, fuel characteristics etc. A key element in analyzing this 
regime's effects are the acoustical resonant modes of the combuster. Indeed 
the first approximation of the high frequency combuster modes includes a 
classical acoustical analysis which is used as a good first approximation to an 
analysis which modifies the classical acoustical analysis to account for the 
effects of combustion. Thus high frequency combustion stability is basically the 
purview of the combuster designer rather than the test stand designer. 
However, to allow a monitoring capability a late version of Dr. Mitchell's (of 
Colorado State sponsored by the USAF under contract F046 1 1 -86-k-0020) 
analysis software has been modified for and installed on the MSFC computers 
in EP-64. Sadly the current state of the art is such that "from scratch" predictive 
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stability analysis is impossible. This is because the detailed 
aerothermochemistry processes involved are not at all well understood. Even 
the equations that are available are complex and difficult to solve. Perhaps with 
time and the renewed interest in liquid propulsion in the USA (after an hiatus of 
twenty years) progress will be made. Various supercomputers e.g. the Cray 
YMP and the TMC Connection Machine (CM) are being employed to investigate 
various single aspects of the aerothermochemical processes. However, much 
is yet to be done even approach simulating an entire combuster with all its 
interconnected processes. 

The state of the art for predicting the lower frequency oscillatory modes is farther 
advanced than that for the high frequency regime. While still depending on a 
number of factors e.g. n and i that define the combustion effects in a from- 
outside-the-black-box or phenomenalistic viewpoint the ability to analyze and 
predict stability questions is reasonably well in hand. And is covered below. 
Unlike the high frequency oscillatory behavior which is almost exclusively the 
province of the combuster designer the lower frequency oscillating modes 
involve intimate coupling between the combuster and the propellant feed 
systems. Thus the lower frequency stability questions are very much the 
concern of the test stand designer (who must of course confer with the 
combuster designers or engine suppliers for certain engine or combuster 
parameters). 

Thus it is seen that except for a much wider use of much more powerful 
computing capacity the state of the art in stability analysis stands not much 
advanced from where it was in 1970. Such manner as those of Von Karman, 
Crocco, Reardon and Tsien are associated with that era (along with many 
others of course, see SP-194). The authors of this report have talked to Dr. 
Reardon during the winter of 1989-90. He had recently set about seeing if a 
sequel to SP-194 could or should be written. He contacted various groups 
once again active in the field and concluded that not enough new technology 
had been developed at this point in time to warrant a new volume on the 
subject. 

Diagrammatically the frequency separation of the modes and stability 
associated phenomena are shown in a general way in the figure below. 
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High Frequency Dynamic 
Processes Confined to 
Combustor 


Diagram Depicting Where Dynamic Processes Occur and Their 

Interaction 
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LOW FREQUENCY STABILITY 


Liquid propellant rocket systems couple the propellant supply system (feed 
ducts, tanks etc) to the combuster through the injection process. Typically a 
number of nozzles are fed from a manifold with fluid pressurized to an average 
pressure substantially above that in the combuster so as to feed the propellant 
into the ignition process at a desired average rate. The specific requirements on 
the injectors and the conditioning of propellant for combustion vary from 
propellant to propellant but for bipropellant non-hypergolic types they may be 
summarized as follows. 

(1) To obtain good atomization of the fluid it is necessary to assure a large 
surface area of contact between the liquid and the hot recirculating gases. 

(2) A proper proportion between the mass of the propellant and the mass of the 
hot gases surrounding it must be obtained. 

(3) A good renewal of the hot gases to activate surface exchanges must be 
assured. 

(4) A fast, uniform as possible, mixing between the two propellants must be 
ensured. 

Some contemplation of the injection and subsequent combustion process 
shows that as the propellants travel downstream in the combustion zone of the 
engine they initially are little effected by the hot gases surrounding them. 
However as they continue down the chamber they begin to mix, be heated by 
the environment and finally begin to burn. Thus in this Crocco style scenario it 
is seen that there is a time delay between the time of propellant injection and 
the time at which it burns or turns from a liquid phase into a gaseous one. 

Crocco et al proposed a model, in wide use, in which the process described 
above is modeled approximately by a time delay or transport lag phenomena 
modeling the time it takes for the liquid propellant to enter the combuster travel 
downstream and (instantaneously) turn into gas. This phenomenon is shown in 
the figure (note the step function approximation). 
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If the time delay from propellant injection to ignition is dominated by one of the 
propellants in a bipropellant system as is the case with, for example, LOX-RP1 
(where the rapid evaporation rate of the LOX causes the RP1 evaporation rate 
to control the combustion process rates) then a single delay model is justified. If 
the delays of both constituents are comparable then a two delay model is to be 
used. Fortunately while it introduces some more computational complexity the 
inclusion of another delay is an easy, straight forward analytical thing to do, as 
will be seen below. 

It is now necessary to consider how variations in combustion parameters 
couples in to the combustion process. This is necessary to establish a feedback 
mechanism or gain between perturbations in the process and the energy 
release due to combustion that can feed an instability if the proper conditions 
occur. Crocco does this in the following way. He assumes that the rate of 
processes at a given location and time are a function of physical factors such as 
pressure, temperature and other factors. He then assumes that all these factors 
influence the pressure and develops a functional relationship as follows. Using 
a Taylor's Series Expansion where over bars denote average values and 
primes denote small perturbations he writes 

f(p,T,z,...) = f(p,T,z) + p|| + tJ| + z§ + ... 

where the partials are evaluated at the respective average value of P» T >--- 
Assuming that T and z are correlated to p then 

T = T(p) and 2 = (p) 
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thus factoring the expression above yields 


« _ ^ r _ = - . 1 73f 3T 3f dz df 

f(p,T,z,...) = f(p,T,z,...|l + P^ ^ 



where the barrel quantity is to be evaluated at the average values for the 
quantities involved. At this point the quantity n is defined as the interaction 


index given by 

n = 

p[3f 3T3f dzdf \ 

[fla P a P aT ^ a P a z ~JJ 


Thus the instantaneous processes rate is written 

f(p,T,Z...)=f(p7 T,Z...)[ 1 +nE] 

P 


If one assumes that the process rate is only proportional to a power of pressure 
then 


foe pn 


from which as before 
f(p) ~ f(p) + P'|^ 


f(p) 

f(p) 

f(p) 


f(p) + p'np 1 ** 1 
oc p 11 + p’np 11 ' 1 

1 + 2 ^ 


p J 


Thus the same result as before results under the assumption of only pressure 
dependence. Therefore, under the assumption that the physical factors are 
correlated one can disregard the explicit effects of all the factors except that of 
the pressure. Crocco and Cheng propose a few examples of physical 
processes which suggest values for n. However historical design values seen 
to be the best source of quantitative data for n's value. 


LOW FREQUENCY STABILITY MODELING 


With the simplified models and concepts described above the tools are in hand 
to develop the first and simplest mode, suitable for stability investigations. The 
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overall approach taken is to work with small perturbations in propellant feed 
rate and combuster pressure. This local linearizations around a nominal 
engine operating point allows the use of all the techniques of linear stability 
theory e.g. Nyquist diagrams and Bode Plots. In what follows primes are used 
to denote perturbation variables and the pressure is normalized (divided by) the 
nominal or average combuster pressure and the mass flow rate by the total 
nominal or average propellant mass flow rate. Both normalizing factors are 
denoted by the average or superbar notation. 

Assuming that both propellants undergo the same time delay C*f) to ignition and 
where b denotes burn, i injection, f fuel, ox oxidizer and t time one writes as 
follows for the propellant mass flow rates. 

riloxb(t) = m ox i(t)(t- Tf) 
rilFb(t) = m Fi (t)(t- Tf) 

The admittance expressions of the feed systems which relate propellant 
delivery rate perturbation to combine for pressure perturbations (as seen from 
the combustion chambers) are developed in explicit from elsewhere in the 
report in terms of feedsystem parameters. Here from the combuster location 
they are expressed as 

o — ~P m ox j 

^OX — — 

m p 

C F = S% 

III P 

or rearranging 

m 0 xi _ -P 
G°x p 
mFi _ ~P 
Gf p 

Equating like quantities 

m'oxi_ _ -P' _ m' Fi 
G ox m P ^Fm 
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These expressions link the propellant feedsystem to the combustion process via 
the chamber measure. 

The use of perturbation analysis has been amply displayed in the foregoing by use of 
the prime notation, Taylor's Series etc. This approach will be enlarged upon below 
but before that the major assumptions governing kM frequency dynamic behavior will 
be listed. 

1 . The gas pressure in the combuster is uniform and oscillates about a mean 
value (i.e. no wave motion is considered until the intermediate frequency 
oscillations are considered). 

2. The temperature of the gas in the chamber is uniform and constant i.e. not 
effected by pressure oscillations. 

3. The time lag from injection of propellant to the burned state is the same for 
all propellant elements. 

4. Propellant feedsystem dynamics including wave motion must be taken into 
account. 

5. Combustion is characterized by a delay. 

As already pointed out the explicit expressions for the feedline admittances are in 
hand. Now the explicit relationships for the combustion process must be derived. 

One relationship follows from the principle of mass conservation. It states that the rate 
of change of the mass of gas in the chamber is equal to the difference between that 
entering (i.e. burned) and that leaving the chamber (i.e. emitted) 

= m b (t) - m^O 

Perturbation quantities are now introduced. 

p = p +p : rrib(t) = nib + m b : rhe = nie + ir^ 

thus 


= : mb(t) = iri b (l +^r) : nie = rhe(l + 5^) 

ot dt m b nie 


now define 


8 ^ 


Substitute into the conservation equations 

T7 dp — '■ ■> "> 

V— = m b + m b - me - m b 
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but in the steady state as much mass on the average is leaving the engine as there is 
entering. Thus 
nib = rife 


and 


V— = mb - 


Normalizing by the average mass flow rate and density values yields 

vPB(p\-m b nj 

m9tl p J m m 


The equation of state for a gas is (R is the gas constant) 
p = pTR 


The perturbation development is 

p = p + p' : p = p + p : T = T + T' 


substitute to obtain (second order infinitesimals are as usual set equal to zero) 

(p + p) = [pT + pT + p T + p T]R 

from which _ _ 

P = pTR : p' =[p T' + Tp] R 

and 

a£ 

dt 

dividing by P produces 

at^p j dt( p j atiT/ 

where the last term on the right hand side of the equation is assumed zero. The 
literature suggests that good results may be obtained from this assumption and that 
appears to be, aside from mathematical tractability, its main justification. 

From elementary theory describing the flow of a compressible gas through a nozzle 
(see reference 8, P. 8-37 equation (47) the mass flow rate of the emitted gas may be 
expressed as 


From this the perturbation analysis proceeds as 
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rfle + me = A 


p + p' 

Vt+t' 

but from perturbation theory 
T «T 

thus from the binomial theorem 

-r- . A I— — T 

mu + rru = D - p-^r + p 

Vt 

from which 

mg T 
4 P 2T 


An alternate ratio in terms of the characteristic rocket velocity (C*) to express the last 
term in the last equation is useful. From the definition of C* 

Q* _ PcAt 
m 


from previous discussion 

Pc 

VT 


m 


thus by substitution 
VT 


C*oc 


let 


C* = C* + C* = b[(t + T')]2 = b[(t 2 



and then 

(T' = if 
C* 2 T 

Using this relationship and the previous development to yield 

* * ■ 

nk = L_£l 
nie P c* 


Making all the substitutions developed to this point yields. 

e^/pLmb' p' t CT' 

^UP / m P c* 


It is now necessary to relate C* to the mixture ratio (r=mo/mf) and chamber pressure to 
continue the development. Evidence in the literature suggests that 
C* = f(p c , r) 
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Using the expression for the total differential produces 

C*’ = | 


tec*) 


tec* 

tepcj 

Pc + 

dr 


however other literature states that the first term on the right is negligible with respect 
to the second. Thus 
ldC* ] 


►' . u 

= r - 


dr I 


Dividing by c* yields the desired ratio 
r* dr Jr 


L 

Now it is necessary to develop an expression for r . The rate at which propellant is 
being burned is given by the sum of the rates at which the oxidizer and the fuel burn. 

m b = m = niob + nifb 


Performing the usual perturbation analysis and making use of the definition of the 
mixture ratio yields (after some tedium). 


r' _(1 +r)m' ob (1 +r) 


m 


m 


™fb 


This expression may be used in conjunction with the expression for the c*’/c* ratio to 
substitute into the differential equation developed earlier to produce the defining 
equation for low frequency stability investigations. 



H{‘ 


, (1+fl 


C* 


tec* 

1 dr , 


mob + 


i(l +T)fdC* 
dr 


1 -^r 




A common way to investigate low frequency stability is by means of a Nyquist plot. To 
transform the equation above into a form suitable for such an investigation recall that 

m’oxb = m 0 i(t -T f ) 
m'fb = rhf,(t -T f ) 


Making use of the LaPlace transform to move to the complex s-plane and making all 
the substitutions indicated above yields the stability equation in a form suitable for the 
Nyquist plot. 


[e c s + 1]’ 


c* 


dr 


ni+f) 

tec*y 

C* 

dr 1 


\ - v 


G f =-l 


which is of the well known form (let s=j w) 
K Qw) = -1 


It is evident that not only does Kflw) encompass a transport lag but also includes 
possibly transcendental functions in the admittance functions (from the duct transfer 
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admittances). However it is well established that this presents no conceptual problem 
to applying the Nyquist procedure (see reference nine, p. 467, section 16). 
Realistically it does, however, require a computer to evaluate the contours of interest. 
One word of caution is in order. When using applications programs care should be 
taken to determine how the terms such as the transport lag or the complex hyperbolic 
functions are evaluated. For example, some programs use a Pade' approximation to 
the complex exponential. This approximation is only good over a limited frequency 
range. Other programs such as the new symbolic programs (e.g. MATHEMATICA) 
evaluate the various exponential functions without amplitude approximation although 
the treatment of phase angles greater than stipulated range(e.g. plus or minus pi) 
should be investigated. 

Inspections of the stability equation suggests a simple extension to the case where 
both the fuel and the oxidizer time lags have to be accounted for separately. Merely 
multiply the transport lag exponential inside the brackets and note that each term is 
delayed. One term applies to the oxidizer, the other to the fuel. Now instead of one 
identical time lag being associated with each of the propellant constituents separate 
time lags may be associated with each constituent. This procedure yields the "double 
time lag model." 

Also note in this development that n (the sensitivity index) does not appear. It would 
appear if the variability of the time delay term is taken into account. Indeed the 
literature suggests that in any but the low frequency investigations it must be taken 
into account. However, it may be neglected in the low frequency regime if the 
variability is small compared to the total time delay (see reference three, section 5.3). 

Some examples of the Nyquist plots to be expected are exhibited in the Appendix 
titled FEEDLINE. In this appendix an increasing number of terms are included in an 
example problem to illustrate their cumulative effect. In this example the feedsystem 
admittance functions were calculated based on a proposed MSFC test stand design. 
Some of the engine parameters e.g.steady state mass flow rates were furnished by 
MSFC. However, other parameters such as the variation of C* with respect to mixture 
ratio were estimated from the general literature. 


INTERMEDIATE FREQUENCY STABILITY MODELING 

As previously indicated the next frequency range of interest is the intermediate one. 
In this regime the effects of combustion coupling with the feedsystem play a part as 
does the variability of the combustion time delay (hence introducing n directly). In 
addition the longitudinal oscillatory ("organ pipe") modes of the combustion chamber 
gases play a significant part in the stability analysis i.e. the chamber pressure etc is 
not considered to be constant throughout the chamber volume as was done in the low 
frequency stability analysis. 

Because intermediate frequency stability in a given combuster propellant feed system 
combination depends upon the two combustion related parameters n and % it is 
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customary to characterize graphically the stability condition of the system in terms of a 
parameter plane rather than a Nyquist plot. This plane (the n versus x plane) contains 
contour separated regions denoting either stable or unstable system operation (in the 
small). Thus proximity to the stability boundary is to be avoided (keeping to the stable 
side, of course). An example of such a parameter plane is shown in the Code 
Development section. 


Included immediately below is a synopsis of the intermediate frequency stability 
investigation and results. An extensive mathematical development leading to these 
results is contained in the Appendix titled "Intermediate Frequency Oscillations in a 
Liquid Propellant Rocket Nozzle." 
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Intermediate Frequency Mode 


Problem Description 

The intermediate frequency mode couples the feed lines, combustion chamber, 
and the nozzle. Lower frequency modes couple only with the feed lines and the 
higher frequency modes couple only with the nozzle. Thus, for the intermediate 
mode, the whole system must be modeled. 

The results of the low frequency mode study are used to determine the 
admittances for the fuel and LOX lines. The flow in combustion chamber is a 
complex phenomena and therefore some simplifying assumptions are made. The 
medium in the chamber consists of reactant gases, liquid oxidizer and fuel, and 
the gaseous products of combustion. This is represented as a two-phase mixture 
comprised of a mass-averaged gas comprising all species and a single mass- 
averaged liquid phase. The flow is inviscid (except for the existence of a 
droplet drag). The gases are thermally and calorically perfect. The liquid 
phase is well dispersed throughout the chamber. The variations of the energy 
(internal plus kinetic) of the liquid are neglected. And, for this analysis, 
there is no heat transfer to the walls. 

The physical setup is illustrated by the schematic: 



The engine consists of an injector face where the liquid propellants enter 
and the throat where the gaseous products exit. The rest of the surfaces are 
solid. Conservation requires that the mass of gaseous products leaving must 
equal the mass of the liquid propellants entering. The momentum of the two must 
also be equal. The energy of the products leaving must be equal to the energy 
of the propellants entering plus the energy due to chemical reaction and change 
of phase. The conservation equations are: 

mass: df/dt + 7*(fV) + S/°l '/bt + U*(/°i.'Vl) = 0 

momentum: d(/°V)/dt + 7-(/°V-V) + 7p = -&(/°l ’VO/dt - 9 - (/°l °Vl -Vl ) 

energy: b(/°e8)/bt + 7*(i°e 8 V) + 5(/°l “ei.s)/dt + 7*(f i/eLsV) + 7*(pV) = E 
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The boundary conditions may be expressed as 

1. For the liquid phase, specify the injection velocity and the 
injection density 

2. The admittance conditions are used as the boundary condition 
at the combustion chamber - nozzle interface 

3. The condition at solid surfaces is the vanishing of the normal 
velocity component 

4. There is no information passed upstream through the throat. 

The mathematical development of the equations and boundary conditions are 
presented in the appendix titled “Intermediate Frequency Oscillations in a 
Liquid Propellant Rocket Nozzle” and will not be presented here. The non- 
dimensional ized equations based on small perturbation theory will be summarized. 

Auxiliary equations: 

Xi = (5 - 1)uuo + (1 + r)(dhi./dr)(M/s)e" st T(Gox - rGf)Poo 
Yi = -upo 

Zi = (1/S)upo + /°t 0 ULO 
Wi = 2uuo 

Mi = M{e-”TW(s) + Ms)]Poo - fKs)po) 

where Ms) = (1 + stt)(Gox + Gf) - effect of total mass flux oscillations 

Ms) = (1 + r)[(r/c*)(5c*/5r) - nrsr](Gox - rGf)/r - effect of mixture 

ratio fluctuations 


PC s) = n(1 - e 8T ) 
po = Poo cosh(sx) 
uo = -( 1/8) Poo sinh(sx) 

fx 


ui = Yi + | { [s(Wi - Xi ) + Mi]cosh[s(x - x’)] 
o 

+ s(Yi + Zi)sinh[s(x - x')]}dx’ 


pi /S = -Wi - 


{ [s(Wi - Xi) + Mi]sinh[s(x - x’)] 


+ s(Yi + Zi)cosh[s(x - x’)]}dx’ 
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Applying the boundary conditions at Lc (as developed in section 3.6 of SP-194) 


u’(Lc) +-A p’(Lc) + € o’ (Lc ) = 0 
u’ = UO + U 1 
p’ = po + pi 


o’ = 00 + 01 = -( 1/u) 


{((5 - 1)/8)po + (1 + r)(dht/dr)(Gox 


- rGr )Pooe _9t T}M exp[-s 


(1/u)dx’ ’ ] dx 1 

X * 


C is generally small compared with A and may be deleted for small perturbations. 
For rockets which have a cylindrical combustion chamber and conical nozzle, A 
may be approximated by M(S - 1)/(2S), where M is the Mach Number at the 
intersection of the chamber and nozzle. 

Parameter plane stability can be investigated by setting A = 0 and regard u 
as an independent variable. Then the boundary equation can be used as a 
relation between two of the engine design parameters or operating parameters 
(e.g. n vs x). Values of u in the intermediate mode may be chosen in the range 

1/tt < u < 1/x. 
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CODE DEVELOPMENT 


Three FORTRAN codes are being delivered. These cover the three major 
modes of oscillation in a liquid propellant rocket: low frequency mode 
(feed line-combustion chamber Interaction), intermediate frequency mode 
(feed line-combustion chamber-nozzle interaction), and high frequency mode 
(combustion chamber-nozzle Interaction). Along with the codes, sample runs 
will also be included on the diskettes. Program listing of each of the codes 
will be furnished. 


Low Frequency Program ( FEEDLINE ] 


The code for the low frequency mode oscillations (FEEDLINE) is based on 
the MATHEMATICA program for the Macintosh as presented in appendix titled 
"A MATHMATICA File for Creating a Model of an Engine Feed System". The 
mathematical derivations are given in this report. The FORTRAN program was 
developed to generate input for the intermediate mode program and to allow the 
use of any PC compatible computer to run the code. Because the FORTRAN code 
is compiled and not interpreted, as with MATHEMATICIA, execution time is 
greatly reduced. The increased speed of the compiled program makes it a 
canidate for the primary code for the low frequency mode. 


The code uses various files for the input and output of data. The input 
may be input by a data file, input interactively, or a combination of both. 

If the input is entered interactively, either completely or as changes to the 
data file, the input information is written to the input file for later runs 
or for use by IMODE, the intermediate frequency mode program. The file names 
used by the program are: 

FUEL.DAT - file with parameters for engine, tank, and fuel 
also with piping configuration, 

LOX.DAT - file with parameters for engine, tank, and lox 
also with piping configuration, 

ADMIT.DAT - output of admittances computed (print file), 

ADMIT. PRN - LOTUS file of admittances computed (plot file). 


The code is interactive with the user, asking questions of the user and 
using the answers to determine what needs to be done. The first question 
asked Is "IS THIS SETUP FOR FUEL OR OXIDIZER? ENTER F OR 0.". This question 
is to determine whether to use FUEL.DAT for fuel or LOX.DAT for oxidizer. 

The other questions are asked of the user to read data from the file, input 
the data interactively, or to modify the data on the file. Once the system is 
completely described, then the admittance calculations are begun. This allows 
a range of frequencies to be run. 


Certain results are output to the screen in order for the user to follow 
the progress of the calculations and possibly request additional frequency 
ranges to be run. Further output is directed to a print file (ADMIT.DAT) and 
a plot file (ADMIT. PRN). The print file may later be directed to a printer 
to obtain a hardcopy. The plot file may be used with LOTUS to obtain plots of 
admittance vs frequency. Either frequency or normalized frequency may be used 
for plotting as both are output to the plot file. A sample plot is included 
below. 
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Intermediate Frequency Program CIMQDE) 

The code for the intermediate mode oscillations (IMODE) is an 
implementation of the equations presented in the appendix titled "Intermediate 
Frequency Oscillations In a Liquid Propellant Rocket Nozzle". The intermediate 
mode couples with the piping (using the results from FEEDLINE) and with the 
nozzle (using the combustion chamber - nozzle interface boundary condition). 
The program was designed to allow the parametric study of two variables at a 
given frequency of oscillation. 

The program allows the user to specify an independent and a dependent 
variable. The dependent variable Is iterated upon until the real part of the 
nozzle boundary condition is within 10 -5 of zero. If it doesn’t converge 
within 10 iterations, the current values are printed out and the user can 
decide whether to continue the iteration. By running through a series of 
values of the Independent variable, data for a parameter plane stability plot 
is created. This may be plotted by using LOTUS. A sample n vs tau curve is 
included below. 
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Five files are used by the program. Two are from FEEDLINE: 

FUEL . DAT - file with parameters for engine, tank, and fuel 

also with piping configuration, 

LOX.DAT - file with parameters for engine, tank, and lox 
also with piping configuration. 

The other three files may have any names. The default names are: 

IMODE.INP - file with parameters for the physical configuration, 
IMODE.OUT - output of calculations (print Tile), 

IMODE.PRN - LOTUS file with parametric values (plot file). 

The user may use the default names or specify them at run time. 

The program interactively allows the user to input values, modify values, 
change independent and dependent variables, change independent variable value, 
or exit program. Also, current values of input parameters may be listed or 
the names and dimensions (e.g. ft/sec) of them may be displayed. 
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The high frequency code was obtained from Dr. Mitchell at Colorado State 
University. The program was written for a VAX computer and has been modified 
to run on a PC. Three modes of oscillation are handled; radial, transverse, 
and longitudinal. No coupling with the piping is considered because of the 
high frequencies the program is designed to analyze. A complete description 
of the program is given in Refs. 5 and 6. 


Several changes were made to the code to get it to run on a PC. In one 
subroutine, several variables were dimensioned by (0:100,0:100). However, in 
the instructions, a maximum of 10 was allowed for the upper dimension. 
Therefore, these were changed to (0:10,0:10) and the modified code would run 
on a PC. Without the changes, the compiled code required more than 640K of 
RAM. Also, a few other aesthetic changes were made. 


Thirteen files are used by the program. The input files are: 
PIN - pressure Input, 

NOZIN - input for nozzle admittance calculation, 

ZONIN - general combustor input, 

CAV1DAT - input for radially oriented absorbers, 

CAV2DAT - Input for axially oriented absorbers, 

RADIN - radial cavity input, 

AXIN - axial cavity input. 

The output files are: 

ZDAT - summary of results and echo of input, 

ZNTAU - summary of stability plot results, 

NTPLOT - contains n, tau data tabulated for plotting, 

RADOUT - results specific to radial absorbers, 

AXOUT - results specific to axial absorbers, 

POUT - contains desired pressure point calculations. 


4 
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Nomenclature 


G - admittance \/f 
a - perturbation in mass flow rate 
P - perturbation in pressure 

( 4 | 

z - impedance 

p - pressure 

R - ratio of pressure drop across a dissipative hydraulic circuit element 

e.g. an orifice to the square of the mass flow rate through the 
element 

r - linearized value of pressure mass flow rate ratio 
_ evaluated at an operating point 

a_ - total propellant flow average value 
- a propellant flow rate average value 
L - hudraulic circuit symbol for inertance 

C - hydraulic circuit symbol for capacitance 

A - cross sectional area of duct 

1 - length of duct 

P - fluid density 

k - bulk modulus of fluid 

PDE - partial differential equation 
S - complex LaPlace operator 
LOX - liquid oxygen 
RP1 - hydrocarbon fuel (kerosene like) 

V - volume of combustion chamber 
T - temperature on an absolute scale 
t - time 

R - gas constant in equation of state 

C - characteristic exhaust velocity 

A - nozzle throat area 

r - mixture ratio (oxidizer mass rate/ fuel mass rate) 
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APPENDIX I 


Four Terminal Network Approximations To Duct Admittances 

In the main body of this report the transfer admittance of a duct is developed 
using a continuous distributed model. This is a good approach especially for 
uniform ducts. Another way to develop a model is to approximate the dynamics 
by assuming a finite element model. This is done by subdividing the duct into a 
number of pieces and then assuming that each piece is characterized by a 
lumped inertance, lumped capacitor and as applicable a lumped resistor. The 
number of elements to choose is an heuristic matter. Aerojet recommends 
choosing an element length eight times smaller than the wavelength of the 
highest harmonic of interest. A more conservative choice would be the oft-used 
order of magnitude. What with powerful computers as ubiquitous as they 
currently are an iterative approach in which successively smaller length 
elements are chosen until some figure of merit e.g. the nth harmonic frequency 
change is met would seem a tractable and rational way to assess the element 
length choice. On notes that in case the duct is not uniform but of arbitrary cross 
section versus position certainly the finite element model technique takes on 
considerable additional appeal. 

To treat the duct or indeed any cascade or tandem connection of elements use 
in made of four-terminal or two terminal pair network theory. This theory is well 
covered in the technical literature e.g. in the works of E.A. Guillemin but enough 
will be developed below for the purposes at hand. Consider the network 
specified in the figure below. 
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The indicated reference conditions (polarities or signs) for pressure and mass 
flow rate are an important part of the basic structure. The network in the box is 
arbitrary, except that it is composed of linear and bilateral elements and does 
not contain any pressure or mass flow rate sources or sinks. In the 


development that follows it is most convenient to use transform notation i.e. The 
LaPlace transform notation which converts time domain dynamics (differential 
equations) to s-domain algebraic ones. A summary of the zero initial condition 
hydraulic element relationship on interest here is given below. 

Capacitor: 



_L 

T 


'w' 

II 

n JL‘ 



P(s ) = Lsm{s) 

-A/V V- 

p(t) = Rm{t) 

P(s) = Rm(s) 


Assuming a topologically flat network in the box it is possible to write the 
following set of equations. 
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Pllttll + P\2 m 2 + + PlnMn = P 1 

P21™1 + P22™2 + - + =-^2 

P3imi + Pl2 m 2 + -. + P3n™n = 0 


The factors Ppq have the dimensions of hydraulic impedance. When p=q they 
are the sum of all impedances on a loop contour. When ,Ppq is + or - the 
impedances Using the "window pane" method of traversing all loops in the 
same direction e.g. clockwise the sign will be negative. As an example 
consider the 3 loop or mesh network 
below 



Pll = Z a + Zfo p\2 = P21 = -Zb P 1 = P a 
P22 = Zb + z c + Zd Pl3 = P31 = 0 P2 = 0 
P33 = Zd + Z e P23 = P32 = -Zd P3 = -Py 

Solving the set of equations consists of determining the k th mass flow rate 
according to 

Mk = Al^ + A^ + + AnkP^ 

A A A 

where k = 1 , 2,... N successively and A is the system determinant 



Pll 

Pl2 

• •• Pin 

A = 

P 21 

P22 

— p2n 


Pn 1 

Pn2 

... Pnn 


and the Ppp terms are the cofactors of A. The actual performance of these 
operations becomes quite tedious if N is greater than about 3. If it should ever 
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be required to determine them for larger order systems a symbolic computer 
language such as MATHEMATICA is suggested. 


Choosing k = 1 and k = 2 yields a solution to the four terminal network problem 
as 


m,=Mi P,-^-P 2 
A A 


\ 


m = Au Pl - !hip 2 
A A 

It is conventional as well as convenient to define the following quantities. 


A n 
yn = — 

yn = yn = 

A 

A22 
y22 = ~ AA 

A 


So that the equations may be written as below. 


m=y n Pi -ynPi 
m 2 = > 21 ^ 1 -y22 p 2 


or 


*il = \y 

?2J 13 ’: 


m 1 
m 


>11 -yn 


21 -y22A lP 2 j 


>1 


These expressions are called the y-system equations, 
from the y coefficients is inverted and left multiplied times 


If the y matrix formed 
the left and right sides 


of the equation above one obtains 

Pi = 2\\m\ - znni2 or [Pil _ \z n -zn] ™i 
P2 = Z21^1 ■ Z22OI2 LP 2J L z 21 -Z22I [m2. 


where 
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z\\ 


yn 

M 


Z\2 = *21 



Z 22 = 


yn 

M 


M = ynyn - yh 


This latter formulation is sometimes referred to as the z-system equations for a 
four terminal network. 


Both the y and z equations are relationships among the four variable 
Pi, Pi, rh\, ni 2 a third set of equations is possible. They may be derived from 
the y and z set as follows. From the second equation in the y set solve for ^1 as 

/>! =Z22p 2 + _L_^ 2 

1 yn 1 yn 1 

Now substitute this expression into the first y equation to yield 


or 


" il =yi + yn ™ 2 ) ' ynP2 


tyl Vi 1 

nil = ^pP 2 + j^rrii ( recall y 12 = ^21) 


A similar procedure carried out on the z set yields 
P 1 = zi±P 2 +tt-m 2 

z 21 Z 21 


• _ In , Z22- 
mi = - — r 2 + - — m 2 

zn Z21 


Because there are now two sets of equations relating the vectors 
[Pi miY and [P 2 mj 7 both can be expressed as the same equation by defining 
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A _ Zjl _y22 
Z 21 >21 

B = ^- = -L- 
Z 21 yn 

c= J^ = Ji 

z 2i y 21 
£> = z 22 _ yi 1 
Z 21 >21 


Thus a simple set of equations may be written in terms of the ABCD parameters 
as 

Pi = AP2 + Bmi 
m\ = CP 2 + Dm2 


It is to be noted that the equations above yield descriptions of the four terminal 
network in terms of three y or three z coefficients whereas the general equation 
has four coefficients. This is because the ABCD are not mutually independent. 
This dependency can be expressed as 


AD - BC = 1 


as can be verified by direct substitution of the y's or z's. 

In as much as the motivation for this Appendix was to treat a duct as a 
succession of finite elements it is now necessary to demonstrate how to 
cascade networks. Actually use of the ABCD methodology renders it 
straightforward. Consider the figure 


For each network the input-output relations are given as 

la-Mta 


Pi 


A' B' 

P '2 

- m l- 


CD'. 

m' 2 - 


however 
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therefore 



= [A P] 

A B' 

P 2 

.mj 1 -C DJ 

.CD'. 

m'x. 


Multiplying out the matrices immediately to the right of the equals mark yields a 
two by two matrix equation relating the input and output pairs of the network. In 
the main body of work input impedance or its reciprocal, admittance, is used. 
Thus an expression for it is required. Given the overall expressions 

Pi = EP 2 + Fm2 
mi = GP 2 + Hrn .2 


The impedance looking into the network at the I - 1' terminals is given by 
_P\ _ EP 2 + Pm2 
nil GP 2 + Hrfi2 


now let 




— Ei. 
m2 


so that 


zi= = 


Ez r +F 
Gz r + H 


is the desired result. EFG and H are obtained from the elements of a matrix 
resulting from multiplying together individual finite element transmission 
matrices. Thus depending upon the desired usage they could either be done 
symbolically or numerically (for a given frequency for example) by computer. 


Perusal of an Aerojet report circa 1967 dealing with some propellant 
feedsystem analyses shows use of yet another alternate formulation. This 
formulation makes use of a propagation function derivable from the work above. 
Closer study shows it to produce a result directly parallel to the development in 
the main body of text in this report. The defining function for the propagation 
function is taken to be 

er=iAD +fBU 
or 

7= InQfAD +mU) 
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Perhaps it is now apparent through the foregoing somewhat forma! 
development that the transfer admittance technique used in the main text is 
highly related to this ABCD technique. As a matter of interest both formalisms 
will be used in the example given below. 

As an example the lumped element model of a duct section element will be 
taken to have the same configuration as that used on an infinitesimal scale for 
the continuous distributed 

model. 


Z2 



Using network reduction techniques as before one writes 


^m = 7}-=Zl , Z2&~ 

Cr2 C/2 

( 7. \/ 7a 1 ^ \ 


22 _ (ziXi + Oiz 2 ) 

Z2 + gt) Cl ( zl + Z2+ cr) 


“ ( z —^) °’( 
G, ( Zi + Z2 + GT> 


G 1 zi(l + G1Z2) 

as the desired transfer admittance ratio for this element. To formulate this 


problem using the ABCD method designate loops and write equations as 
shown below. 
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Z2 



or 


P i = z\m\ - zirii'i 
-P 2 = -zittii + (zi + z 2 )m2 


zi 

. ’Zi 


-z\ 

(Zl + Z2 ) 



These equations are of the form used in the derivation. Therefore by matching 
coefficients 


so 


Zll = Zl Zi2 = Z21 = Zi Z 2 2 — Z\ + Z2 
[z] = Z X Z 2 



B _ [z] _ z l z 2 

Z 21 Z 21 


= z 2 


1 _ 1 D _ Z 22 _ ( Z 1 + Z 2 ) 

Z 21 Z 1 Z 21 Zi 


and from before 


_ Az r + B 
Z ” = Cz r +Z> 


gT +Z2 


1 | ( Z 1 + Z2) 

Z1G1 Zl 


Zl\Z 2 + -^ 


1 

zi + z 2 + 77- 

GiJ 


G 2 = 


zi + z 2 + — 


z l| z 2 + 
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G i z i(l + Z2 G \ ) 

This result is of course precisely the same as previously developed from 
network reduction theory. 

To summarize then if it is desired to cascade finite element models such as 
finite dimensional duct segments modeled as having finite lumped inertance, 
capacitance and resistance several techniques are available. Two are 
developed above. These are the transfer admittance approach as used 
previously and the ABCD approach. The use of either would be greatly 
facilitated in actual practice by use of computer based techniques to perform the 
manipulations (symbolic mathematics) and the numerical evaluations. 
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About FEEDLINE 


The FEEDLINE Mathmatica file can be used to create a model for an 
engine feed system. Using the approach found in document SP -194 by 
Harrje and Reardon, FEEDLINE calculates admittance ratios for each 
component of the feed system. These components may include the 
fuel tank, straight piping, bends in piping, valves, bellows, manifold 
and injector orifices. The file then calculates an overall admittance of 
the system (as seen from the engine), by multiplying these separate 
admittance ratios together. From this equation, frequency response 
plots such as Bode and Nyquist can be created to study low frequency 
instability (chugging), or it can be used as input to the intermediate 
combustion stability model developed separately. 

Shown below is the analytical approach to FEEDLINE, using a simple 
system as in figure 1 . It is important to note where each component 
"occurs" in the feed system, because each admittance ratio depends on 
the admittance ratios of the sections above it. 


The equation for the overall admittance is: 


Go Gi (j 2 G3 G4 G5 G6 G7 Os G9 G10 un 


where. 


Go =admittance at the top of the fuel tank 
Gi 

Go =admittance ratio for the tank 
G2 

Gi =admittance ratio for line 1 


G3 

G 2 =admittance ratio for the valve 

G± 

G 3 =admittance ratio for line 2 
G5 

G4 =admittance ratio for bend 1 
G6 

G 5 =admittance ratio for line 3 
Gj_ 

G6 =admittance ratio for bend 2 
G8 

G7 =admittance ratio for line 4 
G9 

G8 =admittance ratio for the bellows 
G10 

G9 =admittance ratio for line 5 


III 



Gil 

Gio=admittance ratio for the manifold 
Gl2 

Gn=admittance ratio for the orifices 


Since there is no fuel flowing into the top of the tank (for a pressure fed 
system), 

Go=0 

t 

therefore, by cancelation of terms. 


G,=G 1 (pi)(5l)(pi)(pi)(pi)(pl ) (pi)(pi)(5jO)(5u ) (^U ) 
Ul (J2 U3 U4 U5 VJ6 VJ7 CJ8 U9 (jio On 


Each component of the feed system has a general equation for its 
admittance ratio. Therefore, it is easy to accommodate changes in the 
system by simply calculating the additional section’s admittance ratio 
and inserting it in the proper place in the equation above. Adding a 
ratio, however, will cause the subscripts in the following ratios to 
change. Therefore, check all subscripts when an addition is made. 
Another important note to make is that several sections of the feed 
system may have similar parameter values. For example the diameter 
of all the straight lines may be the same, or the bends may have the 
same angles. Therefore, you may want to define a general parameter 
value that may be used in calculating several different admittance 
ratios. Finally, it is very important to use the correct units, or all 
calculations will be inaccurate. The general admittance ratio equations 
for various components, normalized with respect to the combustion 
chamber pressure, are shown below: 


"Admittance ratio for the fuel tank: 


Gjl 

Go=l+sctank 

for Go=0, 

Gi=sctank 

ctank=capacitance associated with the fuel tank 

den*vol*prchamb 
C ktank*tflow 

(% 

den=density of fluid ft 3 
vol=volume of fuel tank (ft 3 ) 

( Jb) 

ktank=bulk modulus of fluid in tank ft 2 
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(ib) 

prchamb=pressure inside combustion chamber ft 2 

/Ibx 

tflow=total flow rate of fluid inside engine ' 
s=jw 

* Assume x&y are subscripts, y=x+l 

* Admittance ratio for constant area feed line: (straight pipe) 

Gy 

G x l+G x Zline tanh(stl) 

Zline=impedance associated with the feed line 

Zline= a*iSsm 

areal*prchamb*g 

a=velodty of sound in the fluid 
areal=area of feed line (ft 2 ) 

(32.2 -ft-) 

g=acceleration of gravity s 2 

tl=time constant of line 

tl=lenl/a (s) 

lenl=length of line (ft) 

^Admittance ratio for orifices: 

G y_ 1 
G x l+ZorG x 

Zor=impedance associated with the orifices 

2*dpror*tflow 

lflow*prchamb 

( ib) 

dpror=pressure drop across orifices ft 2 

dbx 

lflow=flow rate of fluid through the line v s } 

*Note that this method (for orifices) may also be used to find the admittance 
ratio for a valve. 
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* Admittance ratio for bellows: 


model as straight pipe for now 

’‘'Admittance ratio for manifold: (similar to model of fuel tank) 


Gy_ i , sCm 
G x G x 

Cm=capacitance associated with manifold 
(refer to equation used for ctank) 

* Admittance ratio for bends in the line: From Dr. George Doane 

In representing a bend in the feed line, you want to show an 
increase in inertance with constant fluid volume. Therefore, a 
bend is modeled as a straight pipe with a decreased area and an 
increased length. Assume that the bend is two feet long (one 
foot on each side). Therefore, 

11 = 2 (ft) 

A1 = area of straight pipe (ft 2 ) 

Recall, for hydraulic lines: 


L=l/A inertance 
C=(p/k)lA capacitance 


therefore, we can say, 

Ci=#)liA! 

for a constant fluid volume, 
Ci=C 2 

then. 


C 2 =#)l2A 2 

k 

L 2 = a" 

a 2 


liA^bAz 

A 2 


now set L 2 =^Li 

= k. 

where A. is a constant Li 


1 1 - 4 



then. 


a 2 *Aj 

Ail2=XA2ll 

dividing by h: 

Using the information found previously. 


then. 


i 2 =ii^ 

A 2 

A l= X(^2.)A 2 

Ai 

Ai 2 =XA2 2 

Al 

Yx 


a 2 =Ai 


"new area of line" 


Ai "new length of line" 


to find A.: 




where, 



L 2 =L 1 +L" 

L" found from graph in SP194, P. 109 (see APPENDIX A) 
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Using FEEDLINE 


The first step in using FEEDLINE is to break the system up into its 
various components so that each section will have a separate admittance 
ratio. By clearly indicating these divisions, you will avoid confusion in the 
latter calculations. 

One of the most tedious tasks in using this model is obtaining the 
parameter values needed for calculations. One piece of information that is 
important in finding accurate values for viscosity, bulk modulus, etc., is the 
pressure drop down the feed line. Disk 1 is set up to help with these 
calculations. First, however, the analytical approach is discussed below: 

The first step is to calculate the Reynold's number in the line, which 
determines whether the flow is laminar or turbulent. 

Form the equation on p. 81 of Fluid Mechanics : 

jv| R _den*vel*dia 

vis 


(Jt) 

den=density of fluid ft 3 

ft 

vel=velocity of fluid (j) 

v d mow 
den*areal 

where, 

/IK 

lflow=flow rate of fluid in line ' s ' 
areal=area of line (ft 2 ) 


dia=diameter of line (ft) 

(J14 

vis= viscosity of fluid ft' s 

For a Reynold's number greater than 10,000, we can assume that the flow is 
turbulent. 

The next step is to find an approximate friction factor associated with the 
Reynold's number. This can be done using the chart in Appendix B, which 
comes from p. 195 of Fluid Mechanics . Assume that we are using smooth 
pipes. 
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We then want to find the head loss in each section using the equation found 
on p. 93 of Fluid Mechanics : 


head 


_fr*len*vel 2 

dia*2*g 


where, 

fr=friction factor 
len=length of section (ft) 


With this information, we can calculate the pressure drop in the pipe, and the 
pressure drop due to friction. 

For the pressure drop in the pipe we use the equation from p. 98 of Fluid 
Mechanics: 


dpri = d e n*lenl 

v 144 

And, for the pressure drop due to friction, the equation on p. 98 of Fluid 
Mechanics: 


d p rfr= h ead* iien 

144 

With these values for each component of the feed system, we can find the 
total pressure drop, which in turn, allows us to find more accurate values of 
the bulk modulus of the fluid, viscosity of the fluid, etc... 
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On the computer 

Insert DISKI to calculate the pressure drop in the line. Click on to the 
FEEDLINE1 disk icon and then four Mathmatica icons will appear to choose 
from. Each file is set up for various systems. There are two LOX systems, one 
with bends and one without bends, and two RP-1 systems, one with bends 
and one without bends. Choose the appropriate file by clicking onto the icon, 
and Mathmatica will then be started. You are now ready to make changes or 
add new parameters to the file. 

*It is important to remember, when entering these parameters, that you want 
to evaluate each entry as a separate piece of information. To do this, press 
shift/return after you are finished working on the entry, and a bracket should 
appear on the right side of the screen as shown below: 


den = 53.3 


Note that in the example above, the bracket has an extra stem on the upper 
side. This indicates that the cell is INACTIVE. If the entry on the screen 
appears with a bracket containing only one stem rather than two stems as 
shown below. 


den = 53.3 


then this cell is ACTIVE. 

In order for Mathmatica to recognize this cell and use the information in it, 
you want the cell to be ACTIVE. To activate an INACTIVE cell, click onto the 
cell and go up to CELLS. Under CELLS, the word INACTIVE should appear 
with a check mark beside it. Go down to the word INACTIVE and highlight 
it. This will make the check mark disappear, and the cell will be ACTIVE and 
ready to be evaluated by Mathmatica. 

*In all of these files, general parameter abbreviations have been designated 
with meanings and units shown in the list in APPENDIX C. 
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*Once again, it is important that the proper units are used in defining these 
parameters, or the calculations will be inaccurate. 

For the pressure calculations on DISKI, you must first find a value for the 
friction factor. Notice that the file already contains values for the friction 
factors. However, those factors are associated with the parameter values 
initially used, therefore, if any parameter value is changed, you must 
evaluate the notebook again with the new values entered. The easiest way to 
do this is to first make the parameter changes, checking that all those cells are 
ACTIVE. Then, you want to go to the top of the screen under SPECIAL, and 
choose EVALUATE NOTEBOOK. Mathmatica will then evaluate every 
ACTIVE cell in the notebook. You then want to look to see what the new 
Reynold's numbers are. Using these values and assuming that the pipes are 
smooth, use the chart in APPENDIX B to find the associated friction factors. 
Once you have those values, go back to your file and enter them in the proper 
place. Now, go under SPECIAL and pick EVALUATE NOTEBOOK again. 
Mathmatica will evaluate every ACTIVE cell, and will give you the new 
pressure drop calculations. 

^Remember that you must have pressure drop equations for every section of 
the feed system. Make sure that each equation contains the parameters with 
the proper subscripts associated with that section. 

With the information from this file, you can go to various references (see the 
list in APPENDIX D) to find values for the bulk modulus, viscosity, etc... 

Once you have found all the necessary parameters, you are ready to use DISK2 
to calculate the admittance ratios. 

Insert DISK2 and click on to the FEEDLINE2 disk icon. As before, you will see 
four icons, representing the four different systems. Choose the appropriate 
icon and the Mathmatica program will begin. 

You will then see a list of parameters similar to those in the files on DISKI. 

In the same manner as before, you may make changes in these values, 
checking that each cell is still ACTIVE. After making these changes, you 
should once again, check all the equations to make sure that the proper 
subscripts are inserted, and that each section is included. You may then 
evaluate the entire notebook by going up to SPECIAL, and choosing 
EVALUATE NOTEBOOK. Another option is to evaluate each entry step by 
step, by pressing shift/ return for each entry. 

*The Mathmatica file that you are working on may contain plot commands 
from previous runs. Since these plots require a lot of memory and time to 
run, it is best to make these cells INACTIVE when they are not needed, so 
that Mathmatica will ignore them when you choose EVALUATE 
NOTEBOOK. 
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Frequency Plots 


We have found in running these files, the feed systems including the bends 
require a tremendous amount of calculations to produce a plot. This is 
because of the larger number of admittance ratios in these systems. When 
trying to produce frequency plots of these systems, the kernel of Mathmatica 
is usually destroyed due to lack of memory. We realize that the cause of this 
problem stems from the fact that in the initial setup of these files, Mathmatica 
will not remember the previous calculations it has made. Therefore, it 
recalculates every equation over each time a new admittance ratio is figured. 
One can see from this that the time and memory required for calculation as 
admittance ratios are added grows exponentially. Efforts are still being made 
to write commands that will cause Mathmatica to remember the values it has 
already found, thus using less memory and speeding up the time needed. An 
approach that has been used to get frequency plots is discussed next. 

Once you have entered all necessary parameters (using the systems without 
bends), you may plot your results in different forms of frequency response 
plots. These plots, as said before, require a lot of time and memory to run. 
Therefore, it is best that the commands are only evaluated one at a time 
(using shift /return), rather than using EVALUATE NOTEBOOK. 

Now you are ready to make some plots. Some examples are shown. 
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This plot gives you the magnitude of your admittance versus frequency: 


pgi5 - PerametricPlot[{(uJ*tl)/Pi, fibsIgilwllMiu,. 0001 ,1 000}, 
RHesLabel->{“(iu*tl)/Pi Vgi"}] 

gi 



The basic form for a plot of this type is: 

Parametric Plot [{x, y}, {x, x min, x max} ] 
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This plot gives you the admittance ratio in decibels versus w on a log scale: 


pgi2 -ParametricPIotKLogll 0,ti>l ( 20*LogM 0,Rb$[gi[iuIII} > {iu l 62.8,6280}, 


flHesLabel->{ M logliury , db"}] 

db 
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*Note that the previous plots were made using the feed systems without 
bends. Since the time required for plotting the feed systems with bends is so 
great, you may just want to create a table of a few points which can then be 
plotted. An example follows. 

As a note, this table containing only five points required about 30 or 45 
minutes to be produced: 


tgiiub - Table({w*.0202284/3.14,Rbs[gi[iv]]} l {iv, 50, 250,50}] 
{{0.322108, 0.137589}, {0.644217, 0.244106), {0.966325, 0.941706}, 

{1.28843, 0.0210737}, {1.61054, 0.338609}} 
ptgiwb ■ ListPlotltgiwb] 


0 . 8 - 

0 . 6 - 

0.4- 

0 . 2 - 

— i > ■ - t ■ - t * 1 t — 

0.4 0.6 0.8 1.2 1.4 1.6 


-Graphics- 
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Even though the table may not produce a meaningful plot, when compared 
with the plot of the system without bends, you can get a good idea of its form. 
This example uses the show command to put plots on top of each other. 



-Grephics- 
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Nyquist Plots 

Using the overall admittance equations for both the LOX and RP-1 feed 
systems, you may create Nyquist plots to aid in stability checking, etc 

To do this, you first need to get both overall admittance equations into one 
file, in their INPUT format. Insert DISK 2 and get into either the LOX or RP-1 
file. To get an equation, enter, 

gi[w] (shift/ return) 

This will give an equation for the overall admittance. Since the equation will 
initially be in the output format, you may first want to unformat it to save 
space. To do this, click onto the cell and go up to CELLS. Under CELLS, 
highlight FORMATTED and this should make the check mark disappear, and 
the equation will be condensed. Now COPY the equation over. Once it is 
copied, you want the cell to be in the INPUT form by going up to STYLE,going 
down to CELL STYLE and highlighting INPUT. Now, set the equation equal 
to gi in the form: 


gi [wj : = equation 

COPY the equation again in its INPUT form onto the clipboard. Close the 
file on DISK 2, eject DISK 2 and then insert DISK 3. Open up one of the 
Nyquist files and PASTE the equation into the file. You will need to rename 
the equation to something unique, such as 

gilox[w_] or girpl [wj 

since both admittance equations are called gi[w_] on DISK 2. Save the file, 
close it, and eject DISK 3. Insert DISK 2, once again, get into the other feed 
system file, and as before, COPY the overall admittance equation (in its 
INPUT form) onto the clipboard. Eject DISK 2, insert DISK 3, get into the 
same file as before and PASTE the other equation into the file. Once you 
have renamed the second equation, you are ready to form some Nyquist 
Plots. 

*NOTE: There are other methods of transferring the admittance equations to 
the Nyquist files. One may be to COPY one equation into the same file as the 
other. Then you would copy both equations simultaneously onto the 
clipboard for transfer to DISK 3. The choice of method is up to the user. 

The Nyquist plots are basically just a plot of the real versus imaginary plots of 
an equation. Some examples that have been calculated follows: 
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This plot is for the equation: 


X 1 T>-1 I s 


- \vo 't'r 
•Z xg, J 

\ + 0c. 


where. 


© 




are constants 


mu;.* 

gi[ui_] 2*EHpM*uJ , ' , tautI/(1+(thetac*l*iiJ)) 

ln( 12]: m 

pgi ■ ParametricPlot[{Re[gi[u>]],lm[gi[Lu]]},{w, 0.01 ,7000}, 

RKesLabel -> {"re Vim“}l 


im 



-Graphics- 
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This plot is for the equation: 


X2tu3-l ;s 


e J 
l v 



M 


where. 



0 -> ? ) <•* 


A.'T 


are constants 


In [20 ):- 

pgi2 - ParometrlcPlot[{Re[gi[uj]],lm[gi[iii]]},{uj,0.0 1 ,7000}, 
RHesLabel -> {"re","im"}l 
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A final example shows the plot of the equation: 









ln{23]:- 

pgi3a - PoraiuelricPlotI{ReIgi[ip]J,JiT)IgiIiii]]},{iii,0.01,1000}, 
RnesLabel -> {"reVim"}] 

im 
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The process described in the previous pages was used in modeling the feed 
systems shown in APPENDIX E. Since a detailed explanation has already 
been given, a short example will show the calculation of some of the 
parameters for the LOX feed system: 

for LOX: 

line 1: 

lenll = 17.97 (ft) 

diall = 1.04 (ft) = diameter for all lines 

areall = .852 (ft 2 ) = area for all lines 

line 2: 

,lenl2 =8.138 (ft) 

line 3: 

lenl3 =32.51 (ft) 

line 4 

lenl4 =8.652 (ft) 

line 5: 

lenl5 =4.12 (ft) 

75° bend: (assume 1 foot on each side) 

R2 = .52 (ft) (inner radius of bend) 

Rl = 1.56 (ft) (outer radius of bend) 


R2/R1 = .333 

from chart in APPENDIX A: 


-^=.17 

R 1 -R 2 

L „ _ (.17)(1. 56-52) _ 20 g (f[) 

.o52 

L ‘= L ' = ^2 =2 - 35(f,) 

L 2 =L’ +L"= 2.56 (ft) 


X=ft2=2 i 56 = i.09 
Li 2.35 

11 = 2 (ft) 


(assumed length of bend) 
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Ienb75 


Xh 

(1.09)(2) 
2.18 (ft) 


Areab75 = areal /fk 

= .852/vOD? 

= .816 (ft 2 ) 

From AEROJET DATA: 

flowrateLOX = 1 Aowlox = 2264 lb m/s 

flowrateRpi = 1 flowRpi = 848 lb m/s 

tflow = 1 flowLOX + 1 flowRpi = 3112 lb m/s 

^density and viscosity of fluid found in references in APPENDIX D. 
den = 72.13 (lb/ft?) 

vis = 14.32*10"5(lb/ ft s) 

With this information, you can run the file for pressure drop calculations and 
then find remaining parameters. 


11-20 




II- 1 
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Figure 3.2.2c. — Inrrtiaoe £" due to duct curvature. 
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FRICTION FACTOR 
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Parameters for FEEDLINE 


APPENDIX C 


den 

- 

density of fluid (lb/ft^) 

vis 

- 

absolute viscosity of fluid (lb/ft*s) 

areal 

- 

area of line (ft 2 ) 

areab “ 

- 

area of bend (ft 2 ) 

dial 

- 

diameter of line (ft) 

diab 

- 

diameter of bend (ft) 

lenl 

- 

length of line (ft) 

lenb 

- 

length of bend (ft) 

g 

- 

acceleration of gravity (32.2 ft/ s 2 ) 

a 

- 

velocity of sound in fluid (ft/s) 

lflow 

- 

flowrate of fluid through the line (lb/s) 

tflow 

- 

total flowrate of fluid to chamber (lb/s) 

voltank 

- 

volume of tank (ft^) 

volman 

- 

volume of manifold (ft^) 

prchamb 

- 

combustion chamber pressure (lb /ft 2 ) 

dpror 

- 

pressure drop across orifices (lb/ft 2 ) 

veil 

- 

velocity of fluid in line (ft/s) 

velb 

- 

velocity of fluid in bend (ft/s) 

reynolds 

- 

Reynold's number 

frl 

- 

friction factor in line 

frb 

- 

friction factor in bend 

headl 

- 

head loss through line (ft) 

headb 

- 

head loss through bend (ft) 

dprfr 

- 

pressure drop due to friction (lb /ft 2 ) 

dprlin 

- 

pressure drop in line without friction (lb/ ft 2 ) 

kman 


bulk modulus of fluid at manifold pressure 
(lb/ft 2 ) 

ktank 


bulk modulus of fluid at tank pressure 
(lb/ft 2 ) 

cman 

- 

capacitance associated with manifold (sec) 

ctank 

- 

capacitance associated with tank (sec) 

g 

- 

admittance of section 

r 

- 

admittance ratio of each section 

zline 

- 

impedance associated with line 

zbend 

- 

impedance associated with bend 

zor 

- 

impedance associated with orifice 

tl or tb 

- 

time constant (sec) 

gi 

- 

overall admittance seen from engine 


II - IV 



APPENDIX D 


References for FEEDLINE Parameters 

1. AEROJET, "Characterisctics of RP-1 Rocket Fuel," Sacremento, 
California, Libraray Aerojet-General Corporation Liquid Rocket Plant, 
February 14, 1957 

2. Binder, R. C., Fluid Mechaics, New York, Prentice-Hall Inc., 1950, 361 p 

3. Chemical Propulsion Information Agency (CPIA), "Summary of 
Properties," Laurel, Maryland, The Johns Hopkins University Applied 
Physics Laboratory, November, 1987 

4. Harrje, D. T. and Reardon, F. H., "Liquid Propellant Combustion 
Instability," NASA SP-194 (1972) 

5. McCarty, Robert D. and Weber, Lloyd A., "Thermophysical Properties 
of Oxygen from the Freezing Liquid Line to 600 K for Pressures to 5000 PSIA," 
NBS Technical Note 384 (1971) 

6 Von Doehren, Paul J., "Propellant Handbook," Edwards, California, Air 
Force Rocket Propulsion Laboratory, Research and Technology Division, 
January, 1966 

*Other parameters were obtained from NASA data 
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Intermediate Frequency Oscillations 
in a Liquid Propellant Rocket Nozzle 
Wilbur C. Armstrong 


Introduction 


Oscillations In a liquid propulsion rocket nozzle have been of Interest since 
World War II. Possibly they were of Interest before then when liquid rockets 
first began to be studied (e.g. Dr. Goddard). During the 50’s and 60’s much 
progress was made In understanding some of the problems of low frequency and 
intermediate frequency oscillations. In 1941 von Karman suggested Introducing a 
time lag as way to explain combustion Instabilities. Crocco (ref. 1) and others 
developed and applied the idea of a combustion time lag. This use of a time lag 
Is used In the analysis of this report. 

A volume edited by Harrje and Reardon (ref. 2) was published In 1972 which 
presented much of the work done on oscillations in a liquid rocket (SP-194). 
Chapters 3, 4, and 5 form the basis for this paper. An effort will be made to 
develop the fluid dynamic equations and fill In some of the missing mathematics. 
Also, SP-194 has several typographical errors which are difficult to find without 
a detailed following of the mathematics. The emphasis will be on the intermediate 
mode of oscillation; however, the initial development is applicable to all modes. 

In the 1970’s, much work was done by Mitchell and associates at Colorado 
State University with high frequency instability. A report by Mitchell and Eckert 
(ref. 3) describes a simplified computer program for predicting linear 
instabilities. In the early eighties their work slowed down. However, it has 
resumed and Mitchell has a new computer program which should be released soon. 

At an AGARD meeting In October, 1988, Cullck (ref. 4) presented an overview 
paper on combustion Instabilities in liquid propulsion systems. Included Is an 
excellent bibliography for rockets, thrust augmentors, ramjet engines, and passive 
and active control of instabilities. 


m-l 



Problem Description 


The flow In a rocket nozzle Is a complex phenomena and therefore some 
simplifying assumptions will be made. The medium In the chamber consists of 
-eactant gases, liquid oxidizer and fuel, and the gaseous products of combustion. 
This will be represented as a two-phase mixture comprised of a mass-averaged gas 
comprising all species, identified by no subscript, and a single mass-averaged 
liquid phase Identified by the subscript ( )l. The flow is invlscld (except for 
the existence of a droplet drag). The gases are thermally and calorlcally 
perfect. The liquid phase Is well dispersed throughout the chamber. The 
variations of the energy (Internal plus kinetic) of the liquid are neglected. 

And, for this analysis, there Is no heat transfer to the walls. 

The assumptions made allows the use of the following: 

equation of state: p = fRT 


speed of sound a = / 5p/f = / 8RT 

The rocket Is considered to be comprised of two parts: the combustion chamber 
where the Mach number Is low and the nozzle where the Mach number Increases to 1.0 
at the throat. The effect of the fuel lines will be Initially Ignored, but will 
be added later in the discussion. The physical setup is as follows: 


Injector 

face 



The rocket motor consists of an Injector face where the liquid propellants enter 
and the throat where the gaseous products exit. The rest of the surfaces are 
solid. Conservation requires that the mass of gaseous products leaving must equal 
the mass of the liquid propellants entering. The momentum of the two must also 
be equal. The energy of the products leaving must be equal to the energy of the 
propellants entering plus the energy due to chemical reaction and change of phase. 
The conservation equations are: 

mass: df/dt + U-(fV) + 5 /°l Vdt + S-(fi*ViJ = 0 (1) 

momentum: d(/°V)/& t + 9«(/°V-V) + Vp = 'V L )/5t - *M/°l *Vl *Vl) (2) 

energy: b(fe 8 )/d t + 7-(/»e.V) + b(/°L ’ei B )/dt + -(fi. *6 l.V) + (MpV) = E (3) 
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The boundary conditions may be expressed as 

1. For the liquid phase, assign the Injection velocity and the 
Injection density 

2. The admittance conditions are used as the boundary condition 
at the combustion chamber - nozzle Interface 

3. The condition at solid surfaces is the vanishing of the normal 
velocity component 

4. There is no information passed upstream through the throat. 

The last boundary condition is more conveniently used by computing the nozzle 
admittance conditions at the chamber-nozzle Interface and applying the conditions 
there. 


Fluid Flow Equations for Liquid Rocket 

The conservation equations (1-3) are used to study the oscillations In a 
liquid propellant rocket. First the equations are written In a different form. 
Then they are non-dimensional 1 zed. Next, small perturbation theory Is applied and 
finally the one-dimensional longitudinal equations are solved. 

The conservation of mass equation may be rewritten and M defined. 

b/°/bt + V-(/°V) = M = -b/° L ”/bt - V-(/“l°Vl) (4) 

Expand the momentum equation 

V(bf/bt + 7-(/°V)) + fbV/bt + /°V* VV + Vp = 

-Vl (b/’i °/bt - V-(fi/V L )) - Pi °bVc/bt - /“l’Vl-VVl 

and substitute M where appropriate. 

Z°bV/bt + /°V-7V + 7p = M(Vl - V) - Pi 'bVi/bt - /°l'Vl*7Vl (5) 

The next step In the development of these equations Is to non-dimenslonallze 
by use of reference values. Pressure, density, temperature, etc. are divided by 
their values at the Injector face. Velocities are divided by the speed of sound 
at the Injector face, distances by the length of the combustion chamber, and time 
by (length of combustion chamber) / (speed of sound at injector face). Thus 

P = Pdln/pr, U = Udim/ar, X = Xdln/Lc, t - td1»/(Lc/a r), 

M = MdW(/°rar/Lc) , e = edim/hr, etc. 
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Substituting Into equation (4) we obtain 

5(/°/°r)/d(tLc/ar) + S-(/7VVar) = M(/°rar/Lc) = -&(/°L */°r )/5(tLc/a r ) - 7* (Pi YrV L ar ) 
(/°rar/Lc)[d/°/dt + 7-(/°V)] = M(/°rar/Lc) = -(/°rar/Lc)[d/°L °/5t + U-(/°L °Vl)] 

5/°/5t + 7-(/°V) = H = -a^iVat - 9 • (/°l °Vl ) . (6) 

Substituting into equation (5) we obtain 
PPri)( V3r )/5(tLc/3r ) + PPrV&r *7(Var ) + 7*(ppr) 

= M(|»rar/Lc)(VLar-Var) - /°L Vr (5( V L 3r )/&(tLc/ar ) - />L VrV L a r -\7(V L ar ) 
(°rar 2 /Lc{^&V/5t + />V-yV + [pr/(/Var 2 )]9p} 

= />rar 2 /Lc{M(VL-V) - Pi °dV L /dt - /°l“Vl-Wl} 
but, a r 2 = Spr/Z’r giving 

/°5V/5t + /“V-VV + (1/8)7p = M(Vl - V) - m/bt - °Vl -7Vl (7) 

Substituting into equation (3) we obtain 

&(/ ) ^ l reahr )/5(tLc/3r ) + 7 • (^reshrVar )/Lc + &(^l Pr eL s hr )/5(tLc/ar ) 

+ 7*(/°L “/“reLshrVar )/Lc + V • ( PPr Va r )/Lc = EEr 
Pr hr 3r /Lc [b(/*ea )/bt + 7*(/°e 8 V) + &(/°L °e L8 )/5t + 7-(/°L°eL 8 V L )] 

+ (p r ar/Lc )7 • (pV) = EEr 

5(fe.)/bt + V-(pesV) + d(^L°eL 8 )/5t + 7 -(/°l ’eiaVi.) + 7’(pV) = E 

The assumption was made that the internal energy of the liquid did not vary with 
time or space. Also that e 8 = ha - p/P. Thus, 

d(/°e a )/dt + 7-[/*(hs - p/pm + eia[S(Pi 0 )/i>t + 7-(/»l'Vl)] + 7-(pV) = E 
d(fe®)/bt + V-(PhmV) - 7-(pV) - Mslb + 7-(pV) = E 
If e 8 is rewritten as internal energy and chemical energy, then 
5(/°e 8 ) + 7-(/°haV) = MeLa 

The normalization process makes e = T/S, ha = Ta, eLa = 1, T® = T + i (S - DV 2 . 
pea = Pe + i(5 - ^)pv 2 - Pe + p( Ta - T) = PTa + P(e - T) = PJa + /°(T/S- T) 

= pTa + /°T(1/« - 1) = pTa ~ p(8 - 1 )/5 
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Substitute this Into above equation to obtain 

Sin* - ((» - 1 )/8)p]/5t + U*(/»T,V) = H (8) 

Substituting Into the equation of state we obtain 

PPr = PPrRJJr 
PPr = (/°rRTr)/°T 

p = n o) 

Small Perturbation Theory 

If the perturbations of such things as pressure and velocity are small with 
respect to the mean (steady-state) values, the parameters may be represented In 
the form 

p = p + Real(p’e 8t ) (10) 

where p Is the mean value, p’ is the perturbation value, and s is In the form 
s = JS + 1u. The perturbation value Is a function of location only, the time 
variation is taken care of by e 8t . Substituting into equation (6) to obtain 

£>(/“ + /°’e 8 t)/dt + 9-[(f + /°’e 8t )(V + V’e 8t )3 = M + M’e a * 

df/dt + sf’e at + 0-[fV + fV’e at + /°’Ve 8t + / B ’V’e 28t ) = M + Me 8t 

bf/dt + V-(fV) + e^tsf’ + y*(7v’) + y-(f’V) + e 8t \7*(f ’V’ )] = M + M’e 8t 

Subtract H from both sides, divide by e 8t , and set higher order terms to zero. 

S( o’ + U-(fV’) + 9-(/°’V) = M’ (11a) 

similarly, 

— Sf L - 7-(fL*VL’) - \W/*L #, Vl) = M’ (11b) 

Rewrite equation (7) without M and substitute the perturbation equations to 
obtain 

b(fV)/bt + V-iPV 2 ) + (1/8)0*p = -&(fu ’V l )/ bt - 7 • (ft. *V l 2 ) 

b[(f + f’e 8t )(V + V’e 8t )]/5t + V'[(? + /°’e 8t )(V + V’e 8t )(V + V’e 8t )] 

+ (1/8)0-(p + p’e 8t ) = -d[(fL° + Pi ”e 8t )(VL + Vi.’e at )]/dt 
- V-[(pi° + Pi°’e Bt )(yi + Vu’e 8t )(VL + Vi/e 8 *)] 
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d(/°V + /°V’e st + fL’Ve at + /o’ v’e 2st )/&t + 7*(/°V 2 + 2/°W’e 8t + /ov’ 2 e 28t 
+ /°’V 2 e 8t + 2f’ W’e 28t + /o’v’ 2 e 3at ) + (1/8)7*p + (1/8)e 8 t7-p’ 

= ‘Vl ’e 8t + /^‘’VLe 8 * + /°l ° ’Vl ’ e 28t )/dt - V*(/“l°Vl 2 

+ 2/®l ‘VlVl ' e 8t + 'v’L 2 e 2st + /°L a, VL 2 e 8t + 2 /°l ‘’VlVl ’ e 28 t 
+ /°l * ’Vl ’ 2 e 3st ) 

set higher order terms to zero 

(&(/“V) + V • (Tv 2 ) + (1/8)7 *p} + eat[s?V’ + s/°’V + 7*(2/°W’ + f V 2 ) + (1/8)7-p’] 

= - Wl’Vl) + 7*(/°l ° Vl 2 ) } - eat[s^L°VL’ + s/°l°’Vl + V- (2/°l ‘VlVl ’ 

+ * ’Vl 2 )] 

the terms in {} are equal by equation (3) and may be deleted 
sfV* + s/°’v + V-(2fVV’ + /°’V 2 ) + (1/8)7*p’ = -S^l'Vl’ - sf L ” Vl 

- 7* (2/°l ‘VlVl ’ - /°l ° 1 Vl 2 ) 

s(/°V’ + rv) + 7 - (2 /°W’ + /°’V 2 ) + (1/8)V-p’ = -s(/°Z°Vl’ + fL"’Vi) (12) 

- 7* (2/*l ‘VlVl ’ - /°l * ’ Vl 2 ) 

Substitute perturbation equations into equation (8) to obtain 
d[(j» + /°’e 8t )( Ta + Ta ’e 8t ) - ((8 - 1)/8)(p + p’e 8 t)]/5t 

+ 7-[(f + />’e 8t )( Ta + Ta’e 8t )(V + V’e 8t )] = M + H’e 9t 
5[rfa + rra’e 8t + /“’Taeat + /“’Ta’e 28 *:) - ((8 - 1)/8)(p + p’e 8 *)]/5t 

+ 7* [ (^TaV + ^TaV’eat + fTa’Ve 8 * + fTa’V’e 28 * + />’ T 8 Ve 8 t + f’TaV’e 28 * 

+ /*’Ta 1 Ve 2Bt + /a’Ta’V’e 38 *] = M + M’e 8 * 

(Mf Ta - ((8 - 1)/8)p)J/5t + V-(fTaV)} + e 8 t[s(fTa’ + />’Ts - ((8 - 1)/8)p’) 

+ 7 • (^TaV* + fTa’V + /°’TaV)] = M + M’e 8 * 
items in {} are M by equation (8) 

s[?Ta’ + /"’Ts - ((8 - 1 )/8)p’ ] + 7* (^TaV’ + ?Ta ’ V + /° ’ TaV)] = M’ 
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but, T« = 1 and a/’’ + tf-(fV’) + tf*(/°’V) = M’ by equation (11a) 
thus, 

s[?T.’ - ((8 - 1 )/8)p’ ] + (7-(/°VTs’) - 0 (13) 

The energy equation can also be written in terms of entropy (o). Entropy is 
defined by the equation 

p = pd/me-o 

Substituting the perturbation equations gives 
P + p’e at - (p + p’e 8t )( 1/ ^>exp[-(o + o'e Bt )] 

P + = p (i/jf)(i + p’/p e st )( 1/8 )exp(-o)exp(-o’e Bt ) 

Divide by P 

1 + P'/P e Bt = (1 + p’/p e Bt )< 1 /if)exp(-o , e Bt ) 

Replace the factors on the RHS with their series representation and retain only 
up to the first order terms 

1 + p ’ fp = [1 + (1/8)p’/p e 8t ] ( 1 - o’e Bt ) 

1 + P’/P e Bt = 1 - o’e Bt + (1/8)p’/p e Bt - (1/8)p’/po’e 23t ) 

P’/P = -o’ + (1/8)p’/p 

P* p/?= -Pa’ + P 78 

pa’ = p’/8 - TP’ 

Substitute this into equation (13) 

s[/°T’ + (8 - 1)?v.v - ((8 - 1)/8)p’] + V-{MT’ + (8 - 1)V-V’]} = 0 

s[po’ + (8 - 1)?V.V] + 9*{V[po’ + ((8 - 1 )/8)p’ + (8 - DPV-V’D = 0 

s[po’ + (8 - 1)^V-V’3 + 9*{V[p o’ + (8 - 1)?V.V]} = -((8 - 1 )/8)17*(p’V) (14) 

Now apply the perturbation equations to equation (9) 

p + p’e Bt = (P + p’e at )(T + T’e Bt ) 

p + p’e Bt = pf + fT’e Bt + i°’Te Bt + f’T’e 2Bt 

p’ = p’j + pi’ 
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( 15 ) 


which may be rewr-ltten as 
?T* - ((8 - 1)/8)p’ * p’/8 - /*'T 
Also, 

T* + TVe** s T + T’e«t + 4(8 - i)(V + V’e 8t )(V + V’e et ) 

T» + Ts’e 8t s T + T’e 8 * + 4(8 - 1)(W + 2W’e 8t + V’V’e 28t ) 

T» * = T’ + (8 - 1)W’ (16) 

Substitute equation (16) Into equation (13) 

st?T’ + (8 - 1)jW - ((8 - 1)/8)p’) + 5-[i“VT’ + (8 - l)?5*V'] = 0 

use equation (15) to obtain 

s(p75 ~ /“’T) + s(8 - 1 )/ S W*+ 9-[V(/“T’ + (8 - 1)W)] = 0 

s(p’/tf - /°’f) + s(8 - 1)?VV , + 9-[V(p’ - f'T + (8 - 1)fW’)] = 0 

s(r - p’/8) = s( 1 - T)j°’ + s(8 - 1)?W’ + 9- [V(p’ - rr + (5 - 1)?W’)J 

let X = (8 - 1)i®V-V’ + (1 - T )/° ’ 

s(/°’ - p78) = sX + 9- [Vp’ - V/°’T + (8 - DmV.V’)) 

s(7” - p78) = sX + 9-[Vp’ + (1 - T)/*'V + (8 - 1)mv.V’) - /»’V] 

Let Y = -Vp’ + (1 - P)V ’ - (8 - 1)MV-V’) - (1 - T)V/°’ 
s (/>’ - p78) = sX - 9-Y + 9-V’ - 9-(]°V’) - 9 *(A**V) 
or, 

S/°* = Sp’/S + SX - 9-Y + 9-V 1 - 9-(/°V’) - 9-(/°’V) (17) 

Use equation (17) to eliminate f>* from LHS of equations (11a) and (12) 
sp’/S + SX - 9-Y + 9-V’ - 9-(fV’) - 9-(/°’V) + 9-(/>V’) + 9-(/°'V) = M’ 

SP78 + 9-V’ = -sX + 9-Y + M’ (18) 

and now consider equation (12) 

s(?V’ + /°’V) + 9-(27 w’ + f’V 2 ) + (1/8)9-p’ = -s(7l‘Vl’ + fi/’Vt.) 

- 9 -(2/VVi.Vi’ - /°i * ’Vl 2 ) 
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Add sV’ to both sides and rearrange 

sV’ + (1/5)tf-p’ = -s[V/°’ + *Vl ’ + /°l 8, Vl - (1 - f)V’] - 9-[2fW’ + V 2 ] 

- V-(2?l‘VlVl’ - /°l ‘ ’ Vl 2 ) 

Let Z = Vf’ + Fl*Vl’ + Vl/*l * " - (1 - f)V’ 
and W = 2fW’ + 2 ^'VlVl’ + V 2 f’ + Vl 2 ^' 1 
Thus, 

sV’ + (1/5)tf-p’ = -sZ - V-W (19) 

The development of the mass burning rate perturbation is presented In 
section 5.3.1 of SP-194 and will not be given here. The results are 

M’ = H{e-«T[ff(s) + Ms)] - Ms)}p’/p 

where Ms) = (1 + sttKGox + Gf) - effect of total mass flux oscillations 

Ms) = (1 + r)[(r/c*)(5c*/br) - nrSt](Gox - rGf)/r - effect of mixture 

ratio fluctuations 

Ms) = n( 1 - e®*) 

One Dimensional Equations 

When u and u’ are small with respect to the sonic velocity, and thus higher 
order terms may be deleted from the equations. This makes the equations linear In 
terms of the perturbations. Thus, terms such as 

(1 - p), (1 - T), etc. are small and may be deleted. 

This also allows the superposition of solutions. 

Thus, 

p’ = po + pi + 

u’ = uo + ui + 

where all terms of order higher than 1 will be dropped. 
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Consider now the Intermediate mode of oscillation. This mode interacts with 
the feed line system. This mode Is a longitudinal oscillation and we need to 
consider variations In the x - direction only. This simplifies the equations 


zeroth order: (s/8)po + duo/dx = 0 (20a) 

suo + (1/S)dpo/dx = 0 (20b) 

first order: (s/S)pi + dui/dx = -sXi + dYi/dx + Hi (21a) 

sui + (1/8)dpi/dx = -sZi - dWi/dx (21b) 

where Xi * (8 - 1)uuo + (1 + r)(dhi./dr)(H/s)e-* T T(Gox - r0t)Poo (22) 

Yi = -upo (23) 

Zi s (1/b)upo + ^l°uuo (24) 

Wi = 2uuo (25) 

Hi = H{e - » x T[^(s) + Ms)]Poo - tts)po} (26) 


Note: A term has been added to X to account for fluctuating enthalpy (hi.) of the 
liquid propellants resulting from mixture ratio oscillations. 

The zeroth order equations are homogeneous and may be solved 
in a straight forward manner. Use the second equation to remove uo 
from the first equation. This results in an equation in the form 

(D 2 - s 2 )po = 0 

po = A cosh(sx) + B slnh(sx) 

uo = -(1/8) [B cosh(sx) + A slnh(sx)] 

Apply the boundary conditions that at x = 0, p’ = Poo and u’ = 0. 

po = Poo cosh(sx) (27) 

uo = -(1/8)Poo sinh(sx) (28) 

The first order equations are non-homogeneous and are a bit more 
difficult to solve. First obtain an equation In terms of one 
variable. For example, ui . This may be written as 


(D 2 - s 2 )ui = RHS 


where RHS Is the right hand side of the equation. In this case, 
RHS = d(-sXi + dYi/dx + Hi )/dx + s(sZi + dWi/dx) 
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which Is non-homogeneous and may be solved by assuming a solution 
in the form 

ui = vi cosh(sx) + V 2 sinh(sx) 

then take first and second derivatives 

Dui = svi sinh(sx) + sv 2 cosh(sx) 

+cosh(sx)dvi/dx + sinh(sx)dV 2 /dx 

where the second line is set equal to zero. 

D 2 ui = s 2 vi slnh(sx) + s 2 V 2 cosh(sx) 

+s s1nh(sx)dvi/dx + s cosh(sx)dV 2 /dx 

The second line Is equal to RHS because when substituted Into the 
equation 

(D 2 - s 2 )ui = RHS 

everything disappears except the second line. Thus, we have two 
equations in two unknowns. 

cosh(sx)dvi/dx + sinh(sx)dv 2 /dx = 0 

s[s1nh(sx)dvi/dx + cosh(sx)dv 2 /dx] = RHS 

Solving these two equations for dvi/dx and dV 2 /dx 

s[cosh 2 (sx) - s1nh 2 (sx)]dV2/dx = RHS cosh(sx) 

dv 2 = (1/s)RHS cosh(sx)dx 

dvi = -(1/s)RHS s1nh(sx)dx 

thus, 

'X fx 

ui = -cosh(sx)(1/s) RHS s1nh(sx)dx + sinh(sx)( 1/s) RHS cosh(sx)dx 

. O Jo 

Substitute for RHS 

'X 

ui = -cosh(sx)(1/s) [d(-sXi + 7*Yi + Hi )/dx + s 2 Zi + stf*Wi ]s1nh(sx)dx 

. o 

•x 

+ sinh(sx)(1/s) Id(-sXi + \7*Yi + Hi)/dx + s 2 Zi + s7*Wi ]cosh(sx)dx 

. o 
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ui = -cosh(sx)(1/s) 


[d(-sXi + tf*Yi + Mi )/dx]s1nh(sx)dx 


+ s1nh(sx)(1/s) 


[d(-sXi + y*Yi + Mi )/dx]cosh(sx)dx 
o 


- cosh(sx) 


[sZi + 9*Wi ]s1nh(sx)dx + sinh(sx) 
o 


[sZi + 


let u = slnh(sx) and dv = [d(-sXi + tf*Yi + Mi)/dx]dx 


then 


udv = C-sXt + V•y^ + Mt)sinh(sx) - s 
o 


(-sXi + tf-Yi 
o 


substituting Into equation for ui 

ui = -cosh(sx)[(1/s)(-sXi + tf*Yi + Mi)sinh(sx) - 

+ sinh(sx)[(1/s)(-sXi + U-Yi + Mi)cosh(sx) - 


X 

(~SXi 

o 

rx 

(-SXi 

0 


- cosh(sx) 


fx 


[sZi + tf*Wi ]s1nh(sx)dx + slnh(sx) 
o 


[ sZ i + 


ui = cosh(sx) 


(-sXi + 9*Yi + Mi )cosh(sx)dx - slnh(sx) 


rx 

(- 

0 


- cosh(sx) 


[sZi + U*Wi ]s1nh(sx)dx + slnh(sx) 
o 


X 

[sZi + 
o 


17*Wi ]cosh(sx)dx 


+ Mi )cosh(sx)dx 


+ 9*Yi + Mi )cosh(sx)dx] 


+ \7*Yi + Mi )sinh(sx)dx] 


9*Wi ]cosh(sx)dx 


sXi + 7*Yi + Mi )sinh(sx)dx 


V *Wi ]cosh(sx)dx 
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ui = cosh(sx) 


[ s(Wi - Xi ) + Mi + (U*Yi - sWi )]cosh(sx)dx 
o 


- sinh(sx) 


[s(Wi - Xi) + Ht + (9-Yi - sWi)]sinh(sx)dx 
o 


- cosh(sx) 


[s(Yi + Zi ) + (7* Wi - sYi )]s1nh(sx)dx 
o 


- slnh(sx) 


[s(Yi + Zi ) + (U*Wi - sYi )]cosh(sx)dx 
o 


ui = cosh(sx) 


[s(Wi - Xi ) + Mi ]cosh(sx)dx - sinhCsx) 
o 


[s(Wi - Xi) + Mi ]sinh(sx)dx 
o 


- cosh(sx) 


s(Yi + Zi )sinh(sx)dx + sinh(sx) 
o 


s(Yi + Zi )cosh(sx)dx 
o 


+ cosh(sx) 


(\7*Yi - sWi )cosh(sx)dx - sinh(sx) 
o 


(tf*Yi - sWi )sinh(sx)dx 


- cosh(sx) 


(tf‘Wi - sYi )s1nh(sx)dx + slnh(sx) 
o 


(tf*Wi - sYi )cosh(sx)dx 
o 


fx 


ui 


[s(Wi - Xi) + Mi ][cosh(sx)cosh(sx’) - sinh(sx)sinh(sx’ )]dx’ 


s(Yi + Zi )[cosh(sx)s1nh(sx’ ) - sinh(sx)cosh(sx’ )]dx’ 
o 


+ cosh(sx) 


fx 

(9*Yi - sWi )cosh(sx)dx - slnh(sx) 
o 


(9*Yi - sWi )s1nh(sx)dx 
o 


- cosh(sx) 


(\?*Wi - sYi )sinh(sx)dx + slnh(sx) 
o 


(9*Wi - sYi )cosh(sx)dx 
o 
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U1 = 


Us(Wi - Xi ) + Mi]cosh[s(x - x’)l + s(Yi + Zi)s1nh[s(x - x’)]}dx’ + C 
o 


where C = cosh(sx) 


- cosh(sx) 


(U*Yi - sWi )cosh(sx)dx - sinh(sx) 
o 

X 

(9*Wi - sYi )sinh(sx)dx + slnh(sx) 
o 


(\7*Yi - sWi )s1nh(sx)dx 
o 

fx 

(U*Wi - sYi )cosh(sx)dx 
o 


tf*Yi = dYi/dx = -d(upo)/dx = -udpo/dx - podu/dx = -usPoosinh(sx) 

- Poocosh(sx)du/dx 

U*Wt = dWi/dx = 2d(uuo)/dx = 2uduo/dx + uodu/dx = -2u(s/8)Poocosh(sx) 

- (2/8)Poos1nh(sx)du/dx 
substituting into equation for C 


C = cosh(sx) 


[-((8 - 2)/8)sPoou sinh(sx) - Poo cosh(sx)du/dx] cosh(sx)dx 


- sinh(sx) 


- cosh(sx) 


+ slnh(sx) 


[-((8 - 2)/8)sPoou slnh(sx) - Poo cosh(sx)du/dx] sinh(sx)dx 
o 

fx _ _ 

[-((8 - 2)/8)sPoou cosh(sx) - (2/S)Poo s1nh(sx)du/dx] s1nh(sx)dx 


[-((8 - 2)/8)sPoou cosh(sx) - (2/8)Poo s1nh(sx)du/dx] cosh(sx)dx 
o 
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C = -((8 - 2)/8)sPoo{cosh(sx) 


u s1nh(sx)cosh(sx)dx - slnh(sx) 


u s1nh 2 (sx)dx 


+ {cosh(sx) 


u cosh(sx)sinh(sx)dx - sinh(sx) 


-Poo{cosh(sx) 


(du/dx) cosh 2 (sx)dx - slnh(sx) 
o 


u cosh 2 (sx)dx} 


(du/dx)cosh(sx)s1nh(sx)dx 

o 


- cosh(sx)(2/8) 


(du/dx) s1nh 2 (sx)dx + s1nh(sx)(2/8) 
o 


(du/dx)s1nh(sx)cosh(sx)dx) 

o 


C = -((8 - 2)/S)sPoo{cosh(sx) 


u s1nh(2sx)dx - slnh(sx) 


u cosh(2sx)dx} 


-Poo{cosh(sx) 


(du/dx)[1 .+ s1nh 2 (sx)-(2/8)sinh 2 (sx)]dx 


- ((8 - 2)/8) slnh(sx) 


(du/dx) slnh(sx) cosh(sx)dx) 
o 


C = -((8 - 2)/8)sPoo{cosh(sx) 


fx 


u s1nh(2sx)dx - slnh(sx) 


u cosh(2sx)dx) 


-Poocosh(sx) 


(du/dx)dx - ((8 - 2)/8)Poo{cosh(sx) 
o 


rx 

(du/dx)sinh 2 (sx)dx 
o 


- ((8 - 2)/8) slnh(sx) 


(du/dx) slnh(sx) cosh(sx)dx) 
o 


(du/dx)s1nh 2 (sx)dx 

o 


Consider 

let u = s1nh 2 (sx) and dv = (du/dx)dx v = u 


u s1nh 2 (sx) - 2s 


u s1nh(sx)cosh(sx)dx = u sinh 2 (sx) - s 


u sinh(2sx)dx 
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Now, consider 


(du/dx)s1nh(sx)cosh(sx)dx 


let u = s1nh(sx)cosh(sx) and dv = (du/dx)dx v = u 


u s1nh(sx)cosh(sx) - s 


u s1nh(sx)cosh(sx) - s 
thus, 

C = -((8 - 2)/8)sPoo{cosh(sx) 


u sinh 2 (sx)dx - s 


u cosh(2sx)dx 


u cosh 2 (sx)dx 


u sinh(2sx)dx - slnh(sx) 


u cosh(2sx)dx} 


-Poou cosh(sx) - ((5 - 2)/S)Poo{cosh(sx)[u s1nh 2 (sx) 


fx 


- s 

C = Yi 
or, 

ui = Yi + 


u sinh(2sx)dx] - s1nh(sx)[u sinh(sx)cosh(sx) - s 


u cosh(2sx)dx]} 


{ [s(Wi - Xi ) + Mi]cosh[s(x - x’)] 


+ s(Yi + Zi)s1nh[s(x - x’)]}dx’ 
Substituting dui/dx Into equation (21a) 

(s/8) pi + dY i/dx + sWi - sXi + Mi + s{s1nh(sx) 


fs(Wi-Xi) + Mi ]cosh(sx)dx 
o 


- cosh(sx) 


- slnh(sx) 


[s(Wi - Xi ) + Mi ]s1nh(sx)dx + cosh(sx) 


fx 

s(Yi + Zi )cosh(sx)dx 
o 


s(Yi + Zi )s1nh(sx)dx) = -sXi + dYi/dx + Mi 
o 


(29) 
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‘X 


(s/5)pi + sWi = s{ 


[s(Wi - Xi ) + Mi ] [s1nh(sx)cosh(sx’ ) - cosh(sx)s1nh(sx’ )]dx* 


Pi/S = -Wi - 


s(Yi + Zi )[cosh(sx)cosh(sx’ ) - sinh(sx)s1nh(sx’ )]dx’} 
o 


{[s(Wi - Xi ) + Mi]s1nh[s(x - x’)] 
o 


+ s(Yi + Zi)cosh[s(x - x’)]}dx’ (30) 

Applying the boundary conditions at Lc (as developed In section 3.6 of SP-194) 
u’(Lc) + A p’(Lc) + G o’(Lc) = 0 (31) 


where o’ = oo + oi = -(1/u) 


{((5 - 1 )/8)po + (1 + r)(dh L /dr)(Gox 
o 


- rGf )Pooe~ 9T T }M exp[-s 


(1/u)dx”]dx’ 

X ' 


(32) 


C is generally small compared with A and may be deleted for small perturbations. 
For rockets which have a cylindrical combustion chamber and conical nozzle, A 
may be approximated by M(8 - 0/(28), where M is the Mach Number at the 
intersection of the chamber and nozzle. 

When the perturbation equations are substituted into the boundary equation 
(31), a complex equation for s is obtained. The solution for s is difficult. 
However, stability can be Investigated by setting X = 0 and regard u as an 
Independent variable. Then equation (31) can be used as a relation between two of 
the engine design parameters or operating parameters. Values of u In the 
intermediate mode may be chosen between 

1/tt < U < 1/t. 
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Summary 


The equations to analyze Intermediate mode Instabilities have been developed. 
In some cases with long and laborious mathematics. The pertinent equations are 
scattered throughout the report, therefore they have been grouped together for the 


convenience of the reader. 

Xi = (8 - 1)uuo + (1 + r)(dhL/dr)(M/s)e-**T(Gox - rGf)Poo (22) 

Yi = -upo (23) 

Zi = (1/8) upo + 7l°ulo (24) 

Wi = 2uuo (25) 

Hi = Mte-^Tt&s) + Ms)]Poo - fts)po) (26) 

po = Poo cosh(sx) (27) 

uo = -(1/8)Poo sinh(sx) (28) 

*X 

ui = Yi + {[s(Wi - Xi ) + Mi]cosh(s(x - x’)] 

Jo 

+ s(Yi + Zi)s1nh[s(x - x’)]}dx’ (29) 

•x 

pi/8 = -Wi - {[s(Wi - Xi) + Mi]sinh[s(x - x’)] 

. o 

+ s(Yi + Zi)cosh[s(x - x’)]}dx’ (30) 

u’(Lc) + A p’(Lc) + Co’(Lc) = 0 (31) 


o’ = O 0 + 01 = — ( 1/u) 


rx 

{((8 - 1 )/8)po + (1 + r)(dhL/dr)(Gox 
o 


- rGf )Pooe-* T T}M exp[-s 


(1/u)dx”]dx’ 

X ’ 


(32) 


Recall that u’ = uo + ui and p’ = po + pi . 
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Recommendations and Future Plans 

The equations have been developed for obtaining stability limit curves as a 
function of two variables. A computer program may be written for a PC to solve 
the equations for various combinations of variables (e.g. n vs t). If any 
parameter, such as mixture ratio (r), Is changed, new stability limits could be 
computed and compared to the old. 

A FORTRAN program for the PC will be written to determine the stability 
limits as a function of n and t. After the program is checked-out, plans will be 
made to extend the program to allow the use of other parameters to investigate 
stability. 
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Nomenclature 


a - speed of sound 

a r - reference speed of sound 

A - nozzle pressure admittance coefficient 

G - nozzle entropy admittance coefficient 

D - differential operator 

e - static internal energy 

ei.B - stagnation internal energy of liquid 

e 8 - stagnation internal energy e 8 = e + iV 2 dimensional 

e 8 = e + i(8 - DV 2 non-dimensional 

E - energy released per unit time per unit volume in the combustion chamber 
due to chemical reaction and change of phase 

Go x - Oxidizer injection admittance 

Gf - Fuel injection admittance 

hL - enthalpy of liquid 

hs - stagnation enthalpy h 8 = e B + p//° 

Ms) - combustion response function for mixture ratio oscillations 
Ms) - combustion response function for mass flow oscillations 
Lc - length of combustion chamber 

H - combustion response 

Mi - first order term of combustion response perturbation 

n - pressure interaction index 

nr - mixture ratio interaction index 

p - pressure 

pa i in - dimensional pressure 

Pr - reference pressure 
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po - zeroth 'order term of pressure perturbation 
pi - first order term of pressure perturbation 
Poo - maximum pressure at Injection face 
Pis) - combustion response function for pressure 

R - gas constant 

RHS - right hand side of differential equation 
r - mixture ratio 

s - oscillatory term in perturbation equation 

T - temperature 

t - time 

u - axial component of gas velocity 

uo - zeroth order term of velocity perturbation 
ui - first order term of velocity perturbation 
ulo - axial component of liquid velocity 

V - velocity vector of gas 

Vl - velocity vector of liquid 

W - function defined for equation (19) 

Wi - first order term of W 
X - function defined for equation (18) 

Xi - first order term of X 

Y - function defined for equation (18) 

Yi - first order term of Y 

Z - function defined for equation (19) 

Zi - first order term of Z 
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5 - ratio of specific heats 

6 - partial derivative 

* - damping of oscillation 

P - density of gas 

Pl° - mass of liquid per unit chamber volume 

a - entropy 

t - sensitive time lag 

tt - total time lag 

U - frequency of oscillation 
- nab la operator 

Mean values are Indicated by (e.g. p - mean pressure) 

Perturbation values are indicated by ’ (e.g. p’ - perturbation pressure) 
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APPENDIX IV 

FORTRAN Listing of IMODE 
Wilbur C. Armstrong 



Listing of Program IMODE 


PROGRAM MAIN 
C 

C IMODE - INTERMEDIATE MODE OSCILLATIONS 

C 

C VARIABLES IN COMMON /CMPVAL/ COMPLEX 

C 

C XI - FIRST ORDER TERM OF X 

C Y1 - FIRST ORDER TERM OF Y 

C Z1 - FIRST ORDER TERM OF Z 

C W1 - FIRST ORDER TERM OF W 

C Ml - FIRST ORDER TERM OF M 

C PO - ZEROTH ORDER TERM OF PRESSURE 

C PI - FIRST ORDER TERM OF PRESSURE 

C UO - ZEROTH ORDER TERM OF VELOCITY 

C U1 - FIRST ORDER TERM OF VELOCITY 

C RFH - COMBUSTION RESPONSE FUNCTION FOR MIXTURE RATIO 

C RFK - COMPUSTION RESPONSE FUNCTION FOR MASS FLOW 

C RFP - COMBUSTION RESPONSE FUNCTION FOR PRESSURE 


C S - LAMDA + MU i - PERTURBATION OSCILLATION 13 

C GF - FUEL INJECTION ADMITTANCE 14 

C GOX - OXIDIZER INJECTION ADMITTANCE 15 

C RFA - NOZZLE PRESSURE ADMITTANCE COEFFICIENT 16 

C RFC - NOZZLE ENTROPY ADMITTANCE COEFFICIENT 17 

C 
C 

C VARIABLES IN COMMON /RELVAL/ REAL 

C 

C N - PRESSURE INTERACTION INDEX 1 

C TAU - SENSITIVE TIME LAG 2 

C DTAU - DELTA TIME LAG (TOTAL TIME LAG = TAU + DTAU) 3 

C NR - ENTHALPY INTERACTION INDEX 4 

C RBAR - MEAN MIXTURE RATIO 5 

C MBAR - MEAN COMBUSTION RESPONSE FUNCTION 6 

C GAMMA - RATIO OF SPECIFIC HEATS 7 

C POO - MAXIMUM PRESSURE AT INJECTION FACE 8 

C DHLDR - CHANGE IN ENTHALPY WITH CHANGE IN MIXTURE RATIO 9 

C CSTAR - CHARACTERISTIC VELOCITY AT COMBUSTOR EXIT 10 

C DCSDR - CHANGE IN CSTAR WITH CHANGE IN MIXTURE RATIO 11 

C RHOLO - MASS OF LIQUID PER UNIT CHAMBER VOLUME 12 

C ULO - AXIAL COMPONENT OF LIQUID VELOCITY 13 


C LAMDA - DAMPING OF PERTURBATION 
C MU - FREQUENCY OF PERTURBATION 

C TAUT - TOTAL TIME LAG 
C UBAR(50) - VELOCITY ALONG AXIS 
C XBAR(50) - X LOCATIONS ALONG AXIS 

C XLC - X LOCATION OF CHAMBER-NOZZLE INTERFACE (1.0?) 

C 

C 

C VARIABLES IN COMMON /DIMVAL/ REAL 

C 

C ND - PRESSURE INTERACTION INDEX 1 

C TAUD - SENSITIVE TIME LAG SEC 2 

C DTAUD - DELTA TIME LAG (TOTAL TIME LAG = TAU + DTAU) SEC 3 
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NRD - ENTHALPY INTERACTION INDEX 
LAMDAD - OAMPING OF PERTURBATION 
HUD - FREQUENCY OF PERTURBATION 
CDIAM - CHAMBER DIAMETER FT 

TDIAM - THROAT DIAMETER FT 

XLCD - X LOCATION OF CHAMBER-NOZZLE INTERFACE FT 

GAMMAD - RATIO OF SPECIFIC HEATS 

RGAS - GAS CONSTANT FT~2/SEC~2/‘R 

POOD - MAXIMUM PRESSURE AT INJECTION FACE LBF/FT~2 

MBARD - MEAN COMBUSTION RESPONSE FUNCTION LBM/SEC 

RBARD - MEAN MIXTURE RATIO 

DCSDRD - CHANGE IN CSTAR/CHANGE IN MIXTURE RATIO FT/SEC 

DHLDRD - CHANGE IN ENTHALPY/CHANGE IN MIXTURE RATIO FT~2/SEC*2 
RHOLOD - MASS OF LIQUID PER UNIT CHAMBER VOLUME LBM/FT~3 

ULOD - AXIAL COMPONENT OF LIQUID VELOCITY FT/SEC 

PBAR(50) - PRESSURE ALONG AXIS LBF/FT~2 

TBAR(50) - TEMPERATURE ALONG AXIS “R 

XBARD(50) - X LOCATIONS ALONG AXIS FT 


10 

11 

12 

13 

14 

15 

16 

17 

18 


VARIABLES IN COMMON /INTVAL/ 
NVAL - NUMBER OF INPUT POINTS ALONG AXIS 


INTEGER 


PP 

UP 


VARIABLES IN COMMON /RESULTS/ 

P’ = PO + PI 
U’ = UO + U1 


COMPLEX 


c 

SIGP 

- SIG’ = 

SIGO + SIG1 

c 

c 

c 

c 

FUNB 

- BOUNDARY FUNCTION = U’ 

c 

II 

- NUMBER 

OF INDEPENDENT ' 

c 

c 

ID 

- NUMBER 

OF DEPENDENT VAI 

c 


NUMBER 

VARIAB! 

c 


1 

N 

c 


2 

TAU 

c 


3 

DTAU 

c 


4 

NR 

c 


5 

LAM DA 

c 


6 

MU 

c 


7 

CDIAM 

c 


8 

TDIAM 

c 


9 

XLC 

c 


10 

GAMMA 

c 


11 

RGAS 

c 


12 

POO 

c 


13 

MBAR 

c 


14 

RBAR 

c 


15 

DCSDR 

c 


16 

DHL DR 

c 


17 

RHOLO 

c 


18 

ULO 


IV-2 


(o oo n| a (n A 



ooooooooooooooooooooooooooo 


19 PCHMB 

20 TCHMB 

SUBROUTINES 

ADMIT - COMPUTES ADMITTANCE FOR FUEL AND LOX 
BOUND - EVALUATES THE BOUNDARY FUNCTION FUNB 
EVAL - EVALUATES PARAMETERS AT A GIVEN X LOCATION 
FUEL - OBTAINS ADMITTANCE FOR FUEL LINE 

ITER - ITERATES FOR DEPENDENT VARIABLE (REAL BC - ORIGINAL) 
LOX - OBTAINS ADMITTANCE FOR LOX LINE 
NONDIM - NON-DIMENSIONALIZES VARIABLES 
READIN - READS INPUT DATA 
SETVAL - SET VALUES FROM ITERATED VARIABLES 
SETVAR - SET ITERATED VARIABLES FROM VALUES 
ZREAD - READS FREE FORMAT INPUT 

COMPLEX FUNCTIONS 

CCOSH - COMPUTES COMPLEX HYPERBOLIC COSH 
CSINH - COMPUTES COMPLEX HYPERBOLIC SINH 
CTANH - COMPUTES COMPLEX HYPERBOLIC TANH 
FP1 - EVALUATES PI 
FSIGP - EVALUATES SIG’ 

FU1 - EVALUATES U1 


COMMON /CMPVAL/X1 ,Y1 ,Z1 ,W1 ,M1 ,P0,P1 ,U0,U1 ,RFH,RFK,RFP, 

* S , GF , GOX , RFA , RFC 

COMMON /RELVAL/N , TAU , DTAU , NR , RBAR , MBAR , GAMMA , POO , DHLDR , CSTAR , 

* DCSDR , RHOLO , ULO , LAMDA , MU , T AUT , UBAR ( 50 ) , XBAR ( 50 ) , XLC 
COMMON /RESULT/PP,UP,SIGP, FUNB 

COMMON /INTVAL/NVAL 

COMMON /DIMVAL/HOLDD ( 20 ) , XBARD ( 50 ) , PBAR ( 50 ) , TBAR ( 50 ) 

COMMON /TITL/TITLE 

REAL MBAR, N, NR, LAMDA, MU,RVAR(13) 

COMPLEX S,X1 ,Y1 ,Z1 ,W1 ,M1 ,P0,P1 ,U0,U1 , GF , GOX , RFH , RFK , RFP , RFA , RFC 
COMPLEX FP1 , FU1 , FSIGP , PP , UP , SIGP , FUNB , CSINH , CCOSH , CVAR( 1 7 ) 
EQUIVALENCE (N , RVAR( 1 ) ) , (XI , CVAR( 1 ) ) 

CHARACTER*8 VAR(20) , VARP(20) 

CHARACTER* 1 ANS 

CHARACTER* 70 TITLE 

CHARACTER* 2 4 ROCIN,ROCOUT ,ROCPLT 


DATA VAR /’ 

N =’,’ 

TAU =’,' 

DTAU =’,’ 

NR =’,’ 

LAMDA 

* J 

* ’ 

MU =’,’ 

CDIAM =’,’ 

TDIAM =’,’ 

XLC =’,’ 

GAMMA 

* 1 

* ’ 

RGAS =’,’ 

POO = 

MBAR = V 

RBAR 

DCSDR 

“ 1 

* ’ 

DHLDR = 

RHOLO =’,’ 

ULO 

PCHMB = 

TCHMB 

= 7 

DATA VAR P/’ 

N 

TAU 

DTAU ’ , * 

NR *,* 

LAMDA 

> 

f 

* ’ 

MU 

CDIAM ’ , ’ 

TDIAM ’ , ’ 

XLC ’ , ’ 

GAMMA 

» 

> 

* ’ 

RGAS ’ , ’ 

POO 

MBAR ’ , ’ 

RBAR ’ , ’ 

DCSDR 

» 

t 

* i 

DHLDR ’ , ’ 

RHOLO ’ , ’ 

ULO ’ , * 

PCHMB ’ , ’ 

TCHMB 

7 


DATA TOL/. 0001/ 

1 FORMAT (A8 , 1 PE 1 3 . 5 , 2X , A8 , E 1 3 . 5 , ’ FUNB= ’ ,2E13. 5) 

2 FORMAT (A) 

3 FORMAT(/3X,A8,5X,A8, 5X, ’ FUNB(R) ’ ,5X, ’ FUNB(I)V) 

4 FORMAT (1P6E13.5) 
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5 FORMAT ( ’ ’) 

6 FORMAT C ’ “ ’ , A , ’ “ ’ ) 

7 FORMAT C2X, ’ “ * ,A8, ’ “ ’ ,3X, ’ " ’ ,A8, ’ M * ) 

WRITEC*,*) ’ ARE THE FILES YOU ARE USING’ 

WRITEC*,*) ’ IMODE.INP - INPUT DATA’ 

WRITEC*,*) ’ I MODE. OUT - OUTPUT DATA’ 

WRITEC*,*)’ IMODE.PRN - LOTUS PLOT DATA’ 

WRITEC*, ’ (A\) ’ ) ’ ENTER Y OR N ’ 

READ(*,2)ANS 

IF(ANS. EQ. ’Y’.OR.ANS.EQ. ’y’) THEN 
OPEN(5,FILE=’ IMODE.INP’) 

OPEN ( 6 , FILE= ’ IMODE . OUT ’ , STATUS= ’ NEW ’ ) 

OPEN ( 7 , FILE= ’ IMODE . PRN ’ , STATUS= ’ NEW ’ ) 

ELSE 

WRITEC*, ’ (A\) ’ ) ’ ENTER NAME OF FILE CONTAINING INPUT ’ 

READ(*,2)ROCIN 

OPEN(5,FILE=ROCIN) 

WRITEC*, ’ (A\) ’ ) ' ENTER NAME OF FILE FOR OUTPUT ’ 

READ(*,2)ROCOUT 

OPEN ( 6 , FI LE=ROCOUT , STATUS= ’ NEW ’ ) 

WRITEC*, ’(A\)’)’ ENTER NAME OF FILE FOR LOTUS PLOT ’ 
READ(*,2)ROCPLT 

OPEN ( 7 , FI LE=ROCPLT , STATUS= ’ NEW ’ ) 

ENDIF 
XLC=1 .0 
WRITEC*,*)’ ’ 

WRITE(*,*)’ ’ 

WRITEC*,*)’ ’ 

WRITE(*,*) ’ ’ 

WRITE(*,*)’ ’ 

WRITE(*,*)’ ’ 

WRITEC*,*)’ Welcome to IMODE’ 

WRITE(*,*)’ ’ 

WRITE(*,*) ’ Intennediate Mode Rocket Stability Aide’ 

WRITEC*,*)’ ’ 

WRITEC*,*)’ There are three types of input, rocket parameters,’ 
WRITEC*,*)’ Oxidizer feed parameters, and fuel feed parameters,’ 
WRITEC*,*)’ Each may be read from files or from the keyboard’ 
WRITEC*,*)’ ’ 

WRITEC*,*)’ File Name Input’ 

WRITEC*,*)’ ’ 

WRITEC*,*)’ IMODE.INP or NAME read In Rocket Parameters ’ 
WRITEC*,*)’ LOX.DAT Oxidizer Parameters’ 

WRITEC*,*)’ FUEL.DAT Fuel Parameters 

WRITEC*,*)’ ’ 

WRITEC*,*)’ If keyboard entry, you will be prompted for values’ 

GO TO 21 

20 CONTINUE 
WRITEC*,*)’ ’ 

WRITEC*, ’ (A\) ’) ’ Do you want to run another case? Enter Y or N ’ 
READC*,2)ANS 

IFCANS.EQ.’N’.OR.ANS.EQ. ’n’) STOP 

21 CONTINUE 
CALL READIN 
WRITEC*,*)’ ’ 

WRITEC*,*)’ Codes for independent and dependent variables’ 
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WRITE**,*)’ 

CODE 

Variable Name 

CODE 

Variable 

Name’ 

WRITE**,*)’ 

1 

N 

2 

TAU 

1 

WRITE**,*)’ 

3 

DTAU 

4 

NR 

9 

WRITE**,*)’ 

5 

LAMDA 

6 

MU 

9 

WRITE**,*)’ 

7 

CDIAM 

8 

TDIAM 

y 

WRITE**,*)’ 

9 

XLC 

10 

GAMMA 

9 

WRITE**,*)’ 

11 

RGAS 

12 

POO 

9 

WRITE**,*)’ 

13 

MBAR 

14 

RBAR 

9 

WRITE**,*)’ 

15 

DCSDR 

16 

DHLDR 

9 

WRITE**,*)’ 

17 

RHOLO 

18 

ULO 

9 

WRITE**,*)’ 

19 

• 

PCHMB 

20 

TCHMB 

9 


WRITE**,*)’ ’ 

WRITE**, ’(A\)’)’ Enter codes for independent and dependent variabl 
*es ’ 

22 CONTINUE 
READ(*,*)II,ID 
IBAD=0 

IF(II. LT. 1 .OR. II.GT. 20) THEN 
IBAD=1 

WRITE**, *)’ Code for Independent Variable OUT OF RANGE’ 

ENDIF 

IF(I0. LT. 1 .OR. ID.GT. 20) THEN 
IBAD=1 

WRITE**,*)’ Code for Dependent Variable OUT OF RANGE’ 

ENDIF 

IF(IBAD.EQ.I) THEN 
WRITE**, ’ (A\) ’) ’ Enter codes again ’ 

GO TO 22 
ENDIF 

WRITE**,*)’ ’ 

WRITE(6,5) 

WRITE(6,2)TITLE 
WRITE(6,3)VARP(II) ,VARP(ID) 

WRITE(7,6)TITLE 
WRITE(7,7)VARP(II) ,VARP(ID) 

23 CONTINUE 

CALL ITER(ID,TOL) 

WRITE(6, 4)HOLDD(II) , HOLDD(ID) , FUNB 
WRITE(7,4)HOLDD(II) ,HOLDD(ID) 

WRITE**, 1 )VAR(II) ,HOLDD(II) ,VAR(ID) ,HOLDD(ID) ,FUNB 
WRITE**, ’(A\)’) 

* ’ Enter new value for Independent variable (-999 to stop) ’ 

READ* * , * , END=99 ) VAR1 
IF(VAR1.EQ. -999.0) GO TO 20 
CALL SETVAR(VAR1 ,11) 

GO TO 23 
99 CONTINUE 
STOP 
END 
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SUBROUTINE ADMIT (S , GADM , A , AREA , OMAN , CTANK , DPROR , L , LFLOW , PCHMB , 

* SECTN , SEGMN , TFLOW) 

C 

C THIS ROUTINE WILL COMPUTE THE ADMITTANCE RATIO FROM THE TANK TO 

C THE COMBUSTION CHAMBER 

C 

C A = VELOCITY OF SOUND IN THE FLUID (FT/S) 

C AREA(I) = AREA OF FEED LINE (FT2) 

C CMAN = CAPACITANCE ASSOCIATED WITH THE MANIFOLD 

C CTANK = CAPACITANCE ASSOCIATED WITH THE TANK 

C DPROR = PRESSURE DROP ACROSS ORIFICES (LB/FT2) 

C G ( I ) = ADMITTANCE FROM TANK THROUGH PIPE I 

C GADM = OVERALL ADMITTANCE AS SEEN FROM ENGINE 

C GRAV = ACCELERATION OF GRAVITY (32.2 FT/S"2) 

C L ( I ) = LENGTH OF LINE (FT) 

C LFLOW = FLOW RATE OF FLUID THROUGH THE LINE(LB/S) 

C PCHMB = PRESSURE INSIDE COMBUSTION CHAMBER (LB/FT / '2) 

C S = COMPLEX FREQUENCY OF OSCILLATION 

C SECTN(I) = TYPE OF LINE 

C SEGMN = NUMBER OF SECTIONS OF LINE 

C TFLOW = TOTAL FLOW RATE OF FLUID INSIDE ENGINE(LB/S) 

C TL = TIME CONSTANT 

C ZLINE = IMPEDANCE ASSOCIATED WITH THE FEED LINE 

C ZOR = IMPEDANCE ASSOCIATED WITH THE ORIFICES 

C 

COMPLEX CT ANH , G ( 76 ) , GADM , S , W 
REAL AREA(75),L(75), LFLOW 
INTEGER SEGMN, SECTN(75) 

DATA GRAV/32. 2/ 

W=S*A*3. 141593 
G(1)=CTANK*W 
GADM=G(1 )+1 .0 
ZTOP=A*TFLOW/ (GRAV*PCHMB) 

ZOR= 2 . 0*DPR0R*TFL0W/ ( L FLOW*PCHMB ) 

DO 22 I=2,SEGMN+1 
IF(SECTNU-I).EQ.I) THEN 
ZLINE=ZTOP/AREA(I-1) 

TL=L(I-1)/A 

G(I)=(1 .0+CTANH(W*TL)/(G(I-1 )*ZLINE))/(1 .0+G(I-1 )*ZLINE* 

* CTANH(W*TL)) 

GO TO 21 

ELSEIF(SECTN(I-1 ) . EQ. 2) THEN 
G(I)=1.0+(CMAN*W/G(I-1)) 

GO TO 21 

ELSEIF(SECTN(I-1 ) . EQ. 3) THEN 
G(I)=1 .0/(1 .0+Z0R*G(I-1)) 

ENDIF 

21 CONTINUE 

GADM=GADM*G(I) 

G(I)=G(I)*G(I-1 ) 

22 CONTINUE 
RETURN 
END 
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SUBROUTINE BOUND( PP , UP , SIGP , FUNB ) 

COMMON /CMPVAL/X1 ,Y1 ,Z1 ,W1 ,M1 ,PO,P1 ,U0,U1 ,RFH,RFK,RFP, 

* • S , GF , GOX , RFA , RFC 

COMMON /RELVAL/N , TAU , DTAU , NR , RBAR , MBAR , GAMMA , POO , DHLDR , CSTAR , 

* DCSDR , RHOLO , ULO , LAMDA , MU , TAUT , UBAR ( 50 ) , XBAR ( 50 ) , XLC 
COMMON /INTVAL/NVAL 

REAL MBAR, N, NR, LAMDA, MU 

COMPLEX S,X1 ,Y1 ,Z1 ,W1 ,M1 ,P0,P1 ,UO,GF,GOX,U1 ,RFH,RFK,RFP,RFA,RFC 
COMPLEX FP1 , FU1 , FSIGP , PP , UP , SIGP , FUNB , CSINH , CCOSH 
C EVALUATE PP, UP, SIGP, AND FUNB 

P1=FP1 (XLC) 

U1=FU1 (XLC) 

PO=POO*CCOSH(S*XLC) 

U0=- ( 1 . O/GAMMA ) *POO*CSINH( S*XLC ) 

PP=P0+P1 

UP=U0+U1 

SIGP=FSIGP(XLC) 

FUNB=UP+RFA*PP+RFC*SIGP 

RETURN 

END 
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COMPLEX FUNCTION CCOSH(S) 
COMPLEX S 
REAL LAM DA, MU 
LAMQA=REAL(S) 

MU=AIMAG(S) 

COSHR= COSH ( LAMDA ) *COS ( MU ) 

COSHI=SINH(LAMDA)*SIN(MU) 

CCOSH=CMPLX ( COSHR , COSHI ) 

RETURN 

END 


COMPLEX FUNCTION CSINH(S) 
COMPLEX S 
REAL LAMDA, MU 
LAMDA=REAL(S) 

MU=AIMAG(S) 

SINHR=SINH(LAMDA)*COS(MU) 

SINHI=COSH(LAMDA)*SIN(MU) 

CSINH=CMPLX(SINHR , SINHI ) 

RETURN 

END 


COMPLEX FUNCTION CTANH(S) 

COMPLEX S , CTANN , CTAND, CSINH , CCOSH 
CTANN=CSINH(S) 

CTAND=CCOSH(S) 

CTANH=(0. 0,0.0) 

I F ( CTAND . NE . 0 . 0 ) CT ANH=CT ANN/ CT AND 

RETURN 

END 
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SUBROUTINE EVAL(X) 

COMMON /CMPVAL/X1 ,Y1 ,21 ,W1 ,M1 ,P0,P1 ,U0,U1 ,RFH,RFK,RFP, 

* S , GF , GOX , RFA , RFC 

COMMON /RELVAL/N , TAU , DTAU , NR , RBAR , MBAR , GAMMA , POO , DHLDR , CSTAR , 

* DCSDR , RHOLO , ULO , LAMDA , MU , TAUT , UBAR ( 50 ) , XBAR ( 50 ) , XLC 
COMMON /INTVAL/NVAL 

REAL MBAR ,N, NR, LAMDA, MU 

COMPLEX S,X1 ,Y1 ,21 ,W1 ,M1 ,P0,P1 ,U0,U1 ,GF,GOX,RFH,RFK,RFP,RFA,RFC 
COMPLEX CSINH.CCOSH 

C EVALUATE EVERYTHING EXCEPT PP,UP,SIGP 

IF(NVAL.EQ.I) THEN 
UB=UBAR( 1 ) 

GO TO 23 
ENDIF 

DO 21 1=2 ,NVAL 
IF(X. LE.XBAR(I) ) GO TO 22 

21 CONTINUE 
UB=UBAR(NVAL) 

GO TO 23 

22 CONTINUE 

FAC= ( X-XBAR ( I- 1 ) ) / ( XBAR ( I ) -XBAR ( I- 1 ) ) 
UB=UBAR(I-1)+FAC*(UBAR(I)-UBAR(I-1)) 

23 CONTINUE 

RFH=( 1 . O+RBAR ) * ( ( RBAR/CSTAR ) *DCSDR-NR*S*TAU ) * ( GOX 

* -RBAR*GF)/RBAR 

RFK= ( 1 . 0+S*TAUT ) * ( GOX+GF) 

RFP=N*(1.0-CEXP(S*TAU)) 

PO=POO*CCOSH(S*X) 

U0=-( 1 . O/GAMMA )*POO*CSINH(S*X) 

X 1 = ( GAMMA- 1 . 0 ) *UB*UO+ ( 1 . O+RBAR ) *DHLDR* ( MBAR/S ) 

* *CEXP(-S*TAUT ) * ( GOX-RBAR*GF ) * POO 
Y1=-UB*P0 

2 1 = ( 1 . O/GAMMA) *UB*PO+RHOLO*ULO 
W1=2.0*UB*U0 

M 1 =MBAR* ( CEXP ( -S*TAUT ) * (RFK+RFH ) * POO-RFP*PO ) 

RETURN 

END 
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COMPLEX FUNCTION FPI(XL) 

COMMON /CMPVAL/X1 ,Y1 ,Z1 ,W1 ,M1 ,P0,P1 ,UO,U1 ,RFH,RFK, RFP, 

■ * S,GF,GOX,RFA,RFC 

COMMON /RELVAL/N , TAU , DTAU , NR , RBAR , MBAR , GAMMA , POO , DHLDR , CSTAR , 

* OCSDR , RHOLO , ULO , LAM DA , MU , TAUT , UBAR ( 50 ) , XBAR ( 50 ) , XLC 
COMMON /INTVAL/NVAL 

REAL MBAR , N , NR , LAMDA , MU 

COMPLEX S,X1 ,Y1 ,Z1 ,W1 ,M1 ,P0,P1 ,U0,U1 ,GF,GOX,RFH,RFK,RFP,RFA,RFC 
COMPLEX CSINH,CCOSH 
COMPLEX VINT 
C EVALUATE PI 

DX=XL/50.0 
FP1=CMPLX(0. 0,0.0) 

DO 23 1=1,51 
X= ( I— 1 )*DX 
CALL EVAL(X) 

VINT=(S*(W1-X1)+M1)*CSINH(S*(XL-X)) 

* +S*(Y1+Z1)*CCOSH(S*(XL-X)) 

IF(I.EQ.1.0R.I.EQ.51) THEN 

FP1=FP1+0.5»VINT*DX 

ELSE 

FP1=FP1+VINT*DX 

ENDIF 

23 CONTINUE 

FP1=-GAMMA*(W1+FP1 ) 

RETURN 

END 
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COMPLEX FUNCTION FSIGP(XL) 

COMMON /CMPVAL/X1 ,Y1 ,Z1 ,W1 ,M1 ,PO,P1 ,UO,U1 ,RFH,RFK,RFP, 

* S , GF , GOX , RFA , RFC 

COMMON /RELVAL/N , TAU , DTAU , NR , RBAR , MBAR , GAMMA , POO , DHLDR , CSTAR , 

* OCSDR , RHOLO , ULO , LAMDA , MU , T AUT , UBAR ( 50 ) , XBAR ( 50 ) , XLC 
COMMON /INTVAL/NVAL 

REAL MBAR, N, NR, LAMDA, MU 

COMPLEX S,X1 ,Y1 , Z 1 ,W1 ,M1 ,P0,P1 , UO , U1 , GF , GOX , RFH , RFK , RFP , RFA , RFC 
COMPLEX CSINH.CCOSH 
REAL UB(51 ) 

COMPLEX VINT(51 ) , WINT(51 ) , FSIG2 , FCON 

EVALUATE FSIGP (INTEGRATION NOT CHANGED YET) 

DX=XL/50.0 
DO 23 1=1,51 
X=(I-1)*DX 
IF(NVAL.EQ.I) THEN 
UB(I)=UBAR( 1 ) 

GO TO 23 
ENDIF 

DO 21 11=2, NVAL 
IF(X.LE.XBAR(II)) GO TO 22 

21 CONTINUE 

II=NVAL 

22 CONTINUE 

FAC= ( X-XBAR ( I I - 1 ) ) / ( XBAR ( I I ) -XBAR ( II - 1 ) ) 
UB(I)=UBAR(II-1)+FAC*(UBAR(II)-UBAR(II-1)) 

23 CONTINUE 

DO 24 1=1,51 
X=(I-1 )*DX 
CALL EVAL(X) 

V I NT ( I ) = ( ( GAMMA- 1 . 0 ) /GAMMA ) * PO 
VVINT(I)=1.0/UB(I) 

24 CONTINUE 

FCON= ( 1 . O+RBAR ) *DHLDR* ( GOX-RBAR*GF ) *P00 

* *CEXP(-S*TAUT) 

DO 26 1=1,51 

FSIG2=CMPLX(0. 0,0.0) 

DO 25 J=I , 51 

IF(J.EQ.I.OR.J.EQ.51) THEN 
FSIG2=FSIG2+0 . 5*VVINT ( J )*DX 
ELSE 

FSIG2=FSIG2+WINT(J)*DX 

ENDIF 

25 CONTINUE 

FSIG2=CEXP(-S*FSIG2) 

VINT ( I )= (VINT ( I )+FCON)*MBAR*FSIG2 

26 CONTINUE 
FSIGP=CMPLX(0. 0,0.0) 

DO 27 1=1,51 

IF(I.EQ.1.0R.I.EQ.51) THEN 
FSIGP=FSIGP+0. 5*VINT(I )*DX 
ELSE 

FSIGP=FSIGP+VINT(I)*DX 

ENDIF 

27 CONTINUE 
FSIGP=-FSIGP/UB(51 ) 



RETURN 

END 
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SUBROUTINE FUEL(S.GF) 


A 

AREAXI) 

CHAN 

CTANK 

DENS 

DIA(I) 

DPROR 

GF 

KMAN 

KTANK 

L(I) 

LFLOW 

PCHMB 

S 

SECTN(I) 

SEGMN 

TFLOW 

VOL 

VOLMF 


VELOCITY OF SOUND IN THE FLUID (FT/S) 

AREA OF FEED LINE (FT2) 

CAPACITANCE ASSOCIATED WITH THE MANIFOLD 
CAPACITANCE ASSOCIATED WITH THE FUEL TANK 
DENSITY OF FLUID (LBM/FT~3) 

DIAMETER OF LINE (FT) 

PRESSURE DROP ACROSS ORIFICES (LB/FT2) 

OVERALL ADMITTANCE OF FUEL AS SEEN FROM ENGINE 
BULK MODULUS OF FLUID AT MANIFOLD PRESSURE (LBF/FT~2) 
BULK MODULUS OF FLUID IN TANK (LBF/FT~2) 

LENGTH OF LINE (FT) 

FLOW RATE OF FLUID THROUGH THE LINE(LB/S) 

PRESSURE INSIDE COMBUSTION CHAMBER (LB/FT2) 

COMPLEX FREQUENCY OF OSCILLATION 
TYPE OF LINE 

NUMBER OF SECTIONS OF LINE 

TOTAL FLOW RATE OF FLUID INSIDE ENGINE(LB/S) 

VOLUME OF FUEL TANK (FT~3) 

VOLUME OF MANIFOLD (FT~3) 


COMMON /PIPES/PFACE , TFACE 
COMPLEX GF S 

REAL AREA(75),DIA(75),L(75), KMAN, KTANK, LFLOW 
INTEGER SEGMN, SECTN( 75) 

DATA ISTRT/O/ 

1 FORMAT (FI 5.0) 

2 FORMAT ( 15 , 3F1 5. 5) 

IF(ISTRT.EQ.O) THEN 

ISTRT= 1 

OPEN(UNIT= 1 1 , FILE= ’ FUEL . DAT ’ ) 

READ( 11,1) DENS 
READ( 11,1 )TFLOW 
READ(1 1 , 1 )VOLMF 
READ( 11,1) KMAN 
READ( 11,1) PCHMB 
READ (11,1) DPROR 
READ( 11,1 )VOL 
READ( 11,1) LFLOW 
READ( 11,1) KTANK 
READ(11 , 1 )A 
READ(11,1)CTANK 
READ(1 1 , 1 )CMAN 
READ (11, 2) SEGMN 

READ( 11,2) (SECTN(I ) , L( I ) , AREA( I ) , DIA( I ) , 1= 1 , SEGMN) 
ENDIF 

FLOWL=LFLOW*TFACE/TFLOW 

CTANK= ( DENS*VOL*PFACE )/( KTANK*TFACE ) 

CMAN=(DENS*VOLMF*PFACE)/(KMAN*TFACE) 

CALL ADMIT ( S , GF , A , AREA , CMAN , CTANK , DPROR , L , FLOWL , PFACE , 
* SECTN, SEGMN, TFACE) 

RETURN 

END 
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COMPLEX FUNCTION FUI(XL) 

COMMON /CMPVAL/X1 ,Y1 ,Z1 ,W1 ,M1 ,PO, PI ,UO,U1 ,RFH,RFK,RFP, 

* S,GF,GOX,RFA,RFC 

COMMON /RELVAL/N , TAU , DTAU , NR , RBAR , MBAR , GAMMA , POO , DHLDR , CSTAR , 

* DCSDR , RHOLO , ULO , LAMDA , MU , TAUT , UBAR ( 50 ) , XBAR ( 50 ) , XLC 
COMMON /INTVAL/NVAL 

REAL MBAR ,N, NR, LAMDA, MU 

COMPLEX S,X1,Y1,Z1,W1,M1,PO,P1,UO,U1,GF,GOX,RFH,RFK,RFP,RFA,RFC 
COMPLEX CSINH.CCOSH 
COMPLEX VINT 
C EVALUATE U1 

DX=XL/50.0 
FU1=CMPLX(0. 0,0.0) 

DO 23 1=1,51 
X= ( I— 1 )*DX 
CALL EVAL(X) 

VINT=(S*(W1-X1 )+M1 ) *CCOSH (S* ( XL-X ) ) 

* +S*(Y1+Z1 )*CSINH(S*(XL-X) ) 

IFCI.EQ.1.0R.I.EQ.51) THEN 

FU1=FU1+0.5*VINT*DX 

ELSE 

FU1=FU1+VINT»DX 

ENDIF 

23 CONTINUE 
FU1=Y1+FU1 
RETURN 
END 
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SUBROUTINE ITER(ID,TOL) 


ITERATE FOR DEPENDENT VARIABLE 

COMMON /CMPVAL/X 1 , Y 1 , Z 1 , W 1 , Ml , PO , PI , UO , U 1 , RFH , RFK , RFP , 

* S,GF,GOX,RFA, RFC 

COMMON /RELVAL/N , TAU , DTAU , NR , RBAR , MBAR , GAMMA , POO , DHLDR , CSTAR , 

* DCSDR , RHOLO , ULO , LAMDA , MU , TAUT , UBAR( 50) , XBAR( 50 ) , XLC 
COMMON /INTVAL/NVAL 

COMMON /RESULT /PP , UP , SIGP , FUNB 
REAL MBAR , N , NR , LAMDA , MU , RVAR ( 1 3 ) 

COMPLEX S,X1 ,Y1 ,Z1 ,W1 ,M1 ,P0,P1 ,U0,U1 , GF , GOX , RFH , RFK , RFP , RFA , RFC 
COMPLEX FP 1 , FU 1 , FSIGP , PP , UP , SIGP , FUNB , CSINH , CCOSH , CVAR (17) 
EQUIVALENCE (N , RVAR( 1 ) ) , (XI , CVAR( 1 ) ) 

CALL SETVAL(VAL1 , ID) 

CALL BOUND(PP, UP, SIGP, FUNB) 

FUN1=REAL(FUNB) 

IF(ABS(FUN1 ) . LE. TOL) GO TO 22 
VAL2=1.01*VAL1 
IF(VALI.EQ.O) VAL2=0.01 
CALL SETVAR(VAL2, ID) 

CALL BOUND(PP, UP, SIGP, FUNB) 

FUN2=REAL(FUNB) 

IF(FUN1 . EQ. FUN2) THEN 
VAL=VAL1+VAL2 
ELSE 

VAL= VAL 1 - FUN 1 * (VAL2-VAL 1 ) / ( FUN2- FUN 1 ) 

ENDIF 

DO 21 1=1,10 
CALL SETVAR(VAL, ID) 

CALL BOUND(PP, UP, SIGP, FUNB) 

FUN=REAL(FUNB) 

IF(ABS(FUN) . LE.TOL) GO TO 22 
IF(ABS(FUN) . LT.ABS(FUN1 )) THEN 
FUN2=FUN1 
FUN 1= FUN 
VAL2=VAL1 
VAL1=VAL 
ELSE 

FUN1=FUN2 
FUN2=FUN 
VAL 1= VAL 2 
VAL2=VAL 
ENDIF 

IF(FUN1 . EQ. FUN2) THEN 
VAL=VAL1+VAL2 
ELSE 

VAL= VAL 1 -FUN 1 * ( VAL2-VAL 1 ) / ( FUN2- FUN 1 ) 

ENDIF 

21 CONTINUE 

22 CONTINUE 
RETURN 
END 
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SUBROUTINE LOX(S.GOX) 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 


A 

AREA* I) 
OMAN 
CTANK 
DENS 
DIA(I) 
DPROR 
GOX 
KMAN 
KTANK 
L ( I ) 
LFLOW 
PCHMB 
S 

SECTN(I) 

SEGMN 

TFLOW 

VOL 

VOLMF 


VELOCITY OF SOUND IN THE FLUID (FT/S) 

AREA OF FEED LINE (FT"2) 

CAPACITANCE ASSOCIATED WITH THE MANIFOLD 
CAPACITANCE ASSOCIATED WITH THE LOX TANK 
DENSITY OF FLUID (LBM/FT~3) 

DIAMETER OF LINE (FT) 

PRESSURE DROP ACROSS ORIFICES (LB/FT2) 

OVERALL ADMITTANCE OF LOX AS SEEN FROM ENGINE 

BULK MODULUS OF FLUID AT MANIFOLD PRESSURE (LBF/FT2) 

BULK MODULUS OF FLUID IN TANK (LBF/FT^2) 

LENGTH OF LINE (FT) 

FLOW RATE OF FLUID THROUGH THE LINE(LB/S) 

PRESSURE INSIDE COMBUSTION CHAMBER (LB/FT2) 

COMPLEX FREQUENCY OF OSCILLATION 
TYPE OF LINE 

NUMBER OF SECTIONS OF LINE 

TOTAL FLOW RATE OF FLUID INSIDE ENGINE(LB/S) 

VOLUME OF FUEL TANK (FT~3) 

VOLUME OF MANIFOLD (FT~3) 


COMMON /PIPES/PFACE,TFACE 
COMPLEX GOX,S 

REAL AREA(75) , DIA(75) , L(75) , KMAN , KTAN K , L FLOW 
INTEGER SEGMN, SECTN(75) 

DATA ISTRT/0/ 

1 FORMAT ( F 1 5 . 0 ) 

2 FORMAT (I5,3F15.0) 

IF(ISTRT.EQ.O) THEN 

ISTRT=1 

OPEN(UNIT=10, FILE= ’ LOX. DAT’ ) 

READ( 10, 1 )DENS 
READ(10, 1 )TFLOW 
READ(10, 1 )VOLMF 
READ(10, 1 )KMAN 
READ( 10,1) PCHMB 
READ (10,1) DPROR 
READ(10, 1 )VOL 
REAO( 10,1) LFLOW 
READ( 10,1) KTANK 
READ(10, 1 )A 
READ (10,1) CTANK 
READ(10, 1 )CMAN 
READ(10,2)SEGMN 

READ(10,2)(SECTN(I) ,L(I),AREA(I) , DIA(I) , 1=1 , SEGMN) 
ENDIF 

FLOWL=LFLOW*TFACE/TFLOW 

CTANK=(DENS*VOL*PFACE)/(KTANK*TFACE) 

CMAN= ( DENS*VOLMF*PFACE ) / ( KMAN*TFACE) 

CALL ADMIT (S , GOX , A , AREA , CMAN , CTANK , DPROR , L , FLOWL , PFACE , 
* SECTN, SEGMN, TFACE) 

RETURN 

END 
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SUBROUTINE NONDIM(HOLD) 


ROUTINE TO NON-DIMENSIONALIZE INPUT PARAMETERS 

COMtfbN /CMPVAL/X1 ,Y1,Z1,W1,M1,P0,P1 ,U0,U1 f RFH,RFK,RFP, 

* S , GF , GOX , RFA , RFC 

COMMON /RELVAL/N , TAU , DTAU , NR , RBAR , MBAR , GAMMA , POO , DHLDR, CSTAR , 

* DCSDR , RHOLO , ULO , LAM DA , MU , TAUT , UBAR ( 50 ) , XBAR ( 50 ) , XLC 
COMMON /INTVAL/NVAL 

COMMON /DI MVAL/HOL DD ( 2 0 ) , XBARD ( 50 ) , PBAR ( 50 ) , T BAR ( 50 ) 

COMMON /PIPES/PFACE , TFACE 

COMMON /TITL/TITLE 

REAL MBAR , N , NR , LAMDA , MU , RVAR (15) 

REAL MBARD , ND , NRD , LAMDAD , MUD 
REAL HOLD( 20 ) , UBARD( 50 ) , RHOBAR ( 50 ) 

COMPLEX S,X1,Y1,Z1,W1,M1,P0,P1,U0,U1 ,GF,GOX,RFH,RFK,RFP,RFA,RFC 
COMPLEX CVAR(17) 

CHARACTER*8 VAR(13) ,VARD(20) 

CHARACTER* 70 TITLE 

EQUIVALENCE (N , RVAR ( 1 ) ) , (XI , CVAR ( 1 ) ) 

EQUIVALENCE 

* ( ND, HOLDD( 1 ) ) , (TAUD, HOLDD( 2 ) ) , ( DTAUD , HOLDD( 3) ) , 

* (NRD,H0LDD(4)), (LAMDAD, HOLDD(5)), (MUD, HOLDD(6)), 

* (CDIAM,HOLDD(7) ) , (TDIAM,HOLDD(8)) , (XLCD,H0LDD(9) ) , 

* ( GAMMAD , HOLDD ( 1 0 ) ) , ( RGAS ,HOLDD( 1 1 )) , ( POOD ,HOLDD( 12) ) , 

* ( MBARD , HOLDD ( 1 3 ) ) , ( RBARD , HOLDD ( 1 4 ) ) , ( DCSDRD , HOLDD ( 1 5 ) ) , 

* ( DHLDRD , HOLDD ( 1 6 ) ) , ( RHOLOD , HOLDD ( 1 7 ) ) , ( ULOD , HOLDD ( 1 8 ) ) , 

* ( PCHMB , HOLDD( 19) ) , (TCHMB , HOLDD( 20 ) ) 


DATA VAR/’ 

N= ’ , ’ 

TAU= ’ , ’ 

DTAU= ’ , ’ 

NR= ’ , ’ 

RBAR= ’ 

* 

MBAR= ’ , ’ 

GAMMA= ’ , ’ 

P00= ’ , ’ 

DHLDR= ’ , ’ 

CSTAR= ’ 

* 

DCSDR= ’ , ’ 

RHOLO= ’ , ’ 

ULO= V 



DATA VARD/’ 

N = ’,’ 

TAU =’,’ 

DTAU = ’,’ 

NR =’,’ 

LAMDA = 

* 

MU =’,’ 

CDIAM = ’,’ 

TDIAM =’,’ 

XLC =’,’ 

GAMMA = 

* ’ 

RGAS =’,’ 

POO = ’ , ’ 

MBAR = ’ , ’ 

RBAR =’,’ 

DCSDR = 

* 

DHLDR =’,’ 

RHOLO = ’ , ’ 

ULO = ’,’ 

PCHMB =’,’ 

TCHMB = 


DATA PI/3.141 593/ , GC/32 . 1 74/ 

1 FORMAT (A) 

2 FORMAT (A8,1PE13.5,2X,A8,E13.5,2X,A8,E13.5) 

3 FORMAT ( ’ ’) 

N - HOLD(1) 

TAU - HOLD(2) 

DTAU - HOLD(3) 

NR - HOLD(4) 

LAMDA - HOLD(5) 

MU - HOLD(6) 

CDIAM - HOLD(7) 

TDIAM - HOLD(8) 

XLC - HOLD(9) 

GAMMA - HOLD(IO) 

RGAS - HOLD(II) 

POO - HOLD(12) 

MBAR - HOLD(13) 

RBAR - H0LD(14) 

DCSDR - HOLD(15) 

DHLDR - HOLD(16) 
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RHOLO - H0LD(17) 

ULO - H0LD(18) 

PCHMB - HOLD( 19) 

-TCHMB - H0LD(20) 

PBAR - PBAR 
TBAR - TBAR 
XBAR - XBARD 

PCHMB = PBAR ( 1 ) 

TFLOW = LFLOW(LOX) + LFLOW(FUEL) 
LFLOW = LINE FLOW OF LOX OR FUEL 

DO 21 1=1,20 
HOLDD( I )=HOLD( I ) 

21 CONTINUE 

IF(PCHMB.NE.PBAR(1)) THEN 
FAC=PCHMB/PBAR( 1 ) 

DO 22 1=1 , NVAL 
PBAR ( I ) = FAC*PBAR ( I ) 

22 CONTINUE 
ENDIF 

IF(TCHMB.NE.TBARO)) THEN 
FAC=TCHMB/TBAR( 1 ) 

DO 23 1=1 , NVAL 
T BAR ( I ) = FAC*T BAR ( I ) 

23 CONTINUE 
ENDIF 

CAREA=0 . 25*PI*CDIAM**2 
WRITE(6, 3) 

WRITE(6,*)’ CAREA= ’ , CAREA 
TAREA=0. 25*PI*TDIAM**2 
WRITE (6,*) ’ TAREA= ’ ,TAREA 
PFACE=PBAR( 1 ) 

PEXIT= PBAR ( NVAL ) 

TFACE=MBARD 

ASTAR=SQRT(GAMMAD*RGAS*TBAR( 1 ) ) 
WRITE(6,*)’ ASTAR= ’ ,ASTAR 
CSTARD=PEXIT*TAREA*GC/MBARD 
WRITE (6,*)’ CSTARD= ’ ,CSTARD 
DO 24 1=1, NVAL 

RHOBAR ( I )=PBAR ( I ) *GC/ ( RGAS*TBAR (I)) 
WRITE (6,*)’ RHOBAR= ’ , RHOBAR ( I ) 

UBARD ( I )=MBARD/ ( RHOBAR (I)*CAREA) 
WRITE(6,*)’ UBARD=’ ,UBARD(I) 

24 CONTINUE 
N=ND 

T AU=TAUD*AST AR/XLCD 

DTAU= DT AUD* AST AR/XLCD 

TAUT=TAU+DTAU 

NR=NRD 

RBAR=RBARD 

MBAR=MBARD/( RHOBAR ( 1 ) *AST AR*CAREA/XL CD ) 

GAMMA=GAMMAD 

P00=P00D/PBAR( 1 ) 

DHLDR=DHLDRD 
CST AR= CST ARD/ AST AR 
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DCSDR= DCSORD/AST AR 
RHOLO=RHOLOD/RHOBAR ( 1 ) 

ULO=ULGD/ASTAR 
LAHDA= LAMDAD*XLCD/ASTAR 
MU=MUD*XLCD*PI/ASTAR 
XLC=1.0 

DO 25 1=1 ,NVAL 
XBAR(I)=XBARD(I)/XLCD 
UBAR( I )=UBARD( I )/ASTAR 
25 CONTINUE 

S= CHPLX ( LAMDA , MU ) 

CALL FUEL(S,GF) 

CALL LOX(S.GOX) 

RFAR= (GAMMA- 1 . 0)*UBAR( 1 )/(2. 0*GAMMA) 

RFA= CMPLX ( RFAR ,0.0) 

RFC=CMPLX(0. 0,0.0) 

WRITEC*,*) ’ ’ 

WRITEC* , 1 )TITLE 

WRITEC*,*)’ DIMENSIONAL VARIABLES’ 

WRTTFf* NVAIr’’ T51MNVAL 

WRITEC*’, ’(” XBAR= ’ ’ * 1P4E13. 5/(8X,4E13. 5)) ’ )(XBARD(I) , 1=1 ,NVAL) 
WRITEC*, ’ ( ’ ’ UBAR= 1 ’ , 1P4E13. 5/(8X,4E13. 5)) ’ )(UBARD(I) , 1=1 ,NVAL) 
WRITE(*,2)(VARD(I) ,HOLDD(I),I=1 ,20) 

WRITE(6,3) 

WRITE(6, 1 )TITLE 
WRITE(6,3) 

WRITE(6,*) ’ DIMENSIONAL VARIABLES’ 

WRTTF^fi '('* NVAIr’’ T*5V1NVAI 

WRITE (6*, ’ ( ” XBAR= ’ ’ ’, 1P4E13.5/(8X,4E13. 5) ) ’ ) (XBARD(I) ,1=1 , NVAI ) 
WRITE(6, ’ ( ” UBAR= ’ ’ , 1P4E13. 5/(8X,4E13. 5)) ’ ) (UBARD(I) , 1=1 ,NVAL) 
WRITE( 6 , 2 ) ( VARD( I ) , HOLDD( I ) , 1= 1 , 20) 

WRITE(*,*)’ NON-DIMENSIONAL VARIABLES’ 

WRITE(*,’(” NVAL= ” , 15) ’ )NVAL 

WRITEC*, ’(” XBAR= ’ ’ , 1P4E13. 5/(8X,4E13. 5)) ’ ) (XBAR(I) ,1=1 ,NVAL) 
WRITE(*,’(” UBAR= ’ ’ , 1P4E13. 5/(8X,4E13. 5) ) ’ ) (UBAR(I) , 1=1 ,NVAL) 
WRITEC*, ’ ( ” S=” , 1P2E13. 5) ’ ) LAMDA , MU 

WRITE(*,2)(VAR(I) ,RVAR(I) , 1=1,13) 

WRITEC*, ’(” GF= ’ ’ , 1P2E13. 5, 5X, ’ ’ GOX= ” ,2E13. 5) ’ )GF,GOX 

WRITEC*, ’(” RFA=”,1P2E13.5,5X,” RFC= ” ,2E13. 5) ’ )RFA,RFC 

WRITE(6,3) 

WRITE (6,*)’ NON-DIMENSIONAL VARIABLES’ 

WRITE(6, ’ ( ” NVAL= ” , 15) ’ )NVAL 

WRITE(6, ’ ( ” XBAR= ’ ’ , 1P4E13. 5/(8X,4E13. 5) ) ’ ) (XBAR(I) , 1=1 ,NVAL) 
WRITEC6, ’ ( ’ ’ UBAR= ’ ’ , 1P4E13. 5/(8X,4E13. 5)) ’ ) (UBAR(I) , 1=1 ,NVAL) 
WRITE(6, ’ ( ” S= ” , 1P2E13. 5) ’ ) LAMDA, MU 

WRITE(6,2) (VAR(I) ,RVAR(I) , 1=1,13) 

WRITE(6, ’ ( ” GF= ” , 1P2E13. 5,5X, ” GOX= ” ,2E13. 5) ’ )GF,GOX 
WRITE(6, ’ ( ” RFA=”,1P2E13.5,5X,” RFC=” ,2E13.5) ’ )RFA,RFC 

WRITEC*, ’(A\)’)’ Hit ENTER to continue ’ 

READC*,*) 

RETURN 

END 
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SUBROUTINE READIN 


INPUT PARAMETERS 


VARIABLES IN COMMON /DIMVAL/ REAL 


NO - PRESSURE INTERACTION INDEX 1 

TAUD - SENSITIVE TIME LAG SEC 2 

DTAUD - DELTA TIME LAG (TOTAL TIME LAG = TAU + DTAU) SEC 3 

NRD - ENTHALPY INTERACTION INDEX 4 

LAMDAD - DAMPING OF PERTURBATION 5 

MUD - FREQUENCY OF PERTURBATION 6 

CDIAM - CHAMBER DIAMETER FT 7 

TDIAM - THROAT DIAMETER FT 8 

XLCD - X LOCATION OF CHAMBER-NOZZLE INTERFACE FT 9 

GAMMAD - RATIO OF SPECIFIC HEATS 10 

RGAS - GAS CONSTANT FT2/SEC~2/°R 11 

POOD - MAXIMUM PRESSURE AT INJECTION FACE LBF/FT~2 12 

MBARD - MEAN COMBUSTION RESPONSE FUNCTION LBM/SEC 13 

RBARD - MEAN MIXTURE RATIO 14 

DCSDRD - CHANGE IN CSTAR/CHANGE IN MIXTURE RATIO FT/SEC 15 

DHLDRD - CHANGE IN ENTHALPY/CHANGE IN MIXTURE RATIO FT~2/SEC~2 16 
RHOLOD - MASS OF LIQUID PER UNIT CHAMBER VOLUME LBM/FT~3 17 

ULOD - AXIAL COMPONENT OF LIQUID VELOCITY FT/SEC 18 

PBAR(50) - PRESSURE ALONG AXIS LBF/FT~2 

TBAR(50) - TEMPERATURE ALONG AXIS “R 

XBARD(50) - X LOCATIONS ALONG AXIS FT 


* 

* 


* 

* 

* 

* 

* 

* 


* 

* 


COMMON /CMPVAL/X1 ,Y1 ,Z1 ,W1 ,M1 ,P0,P1 ,U0,U1 ,RFH,RFK,RFP, 

S , GF , GOX , RFA , RFC 

COMMON /RELVAL/N , TAU , DTAU , NR , RBAR , MBAR , GAMMA , POO , DHLDR , CSTAR , 
DCSDR , RHOLO , ULO , LAMDA , MU , TAUT , UBAR ( 50 ) , XBAR ( 50 ) , XLC 
COMMON /INTVAL/NVAL 

COMMON /DIMVAL/HOLDD (20) , XBARD( 50 ) , PBAR ( 50 ) , TBAR ( 50 ) 

COMMON /TITL/TITLE 

REAL MBAR , N , NR , LAMDA , MU , RVAR( 1 5 ) 

REAL MBARD , ND , NRD , LAMDAD , MUD , HOLD ( 2 0 ) 

COMPLEX S, XI , Y1, Z1 , W1 ,M1 ,P0, PI ,U0,U1,GF, GOX, RFH,RFK,RFP, RFA, RFC 
COMPLEX CVAR(17) 

EQUIVALENCE ( N , RVAR ( 1 ) ) , ( X 1 , CVAR ( 1 ) ) 

EQUIVALENCE (ND , HOLD( 1 ) ) , (TAUD , HOLD( 2 ) ) , ( DTAUD , HOLD( 3 ) ) , 

(NRD, HOLD(4)), (LAMDAD, HOLD(5)), (MUD, HOLD(6)), 
(CDIAM,HOLD(7)), (TDIAM, HOLD(8)), (XLCD, HOLD(9)), 
(GAMMAD, HOLD( 10)), (RGAS, HOLD(II)), (POOD, HOLD(12)), 
(MBARD, HOLD( 13) ) , (RBARD, HOLD( 14) ) , (DCSDRD, HOLD( 15)) , 
(DHLDRD, HOLD(16)), (RHOLOD, HOLD( 17)), (ULOD, HOLD( 18)), 
(PCHMB,HOLD(19)) , (TCHMB,HOLD(20) ) 

CHARACTERS VAR(20) , VARP(20) ,VARL(20) ,NAME 
CHARACTER* 1 ANS 
CHARACTER*70 TITLE 
DATA IGO/O/ 


ND = V 

TAUD =’ , ’ 

DTAUD = V 

NRD =’ , ’ LAMDAD = ’ 

MUD = ’,’ 

CDIAM = 

TDIAM = ’,’ 

XLCD = V GAMMAD =’ 

RGAS =’,’ 

POOD = V 

MBARD =’,’ 

RBARD =’,’ DCSDRD =’ 


* 


I 
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* ’DHLDRD 

=’,’ RHOLOD 

= ’ , ’ ULOD 

=’,’ PCHMB 

=’,’ TCHMB ='/ 

DATA VARP/’ND 

’,’TAUD 

’ , ’DTAUD 

’ , ’ NRD 

’ , ’ LAMDAD ’ , 

* ’MUD 

’,’CDIAM 

’ , ’TDIAM 

’,’XLCD 

’ , ’ GAMMAD ’, 

* - ’RGAS 

’ , ’POOD 

’,’MBARD 

’,’RBARD 

’,’ DCSDRD ’, 

* ’ DHLDRD 

’, ’RHOLOD 

’ , ’ULOD 

’ , ’ PCHMB 

’, ’TCHMB ’/ 

DATA VARL/’nd 

’ , ’taud 

’ , ’dtaud 

’,’nrd 

’,’lamdad ’, 

* ’mud 

’ , ’cdlam 

’ , ’tdlam 

’,’xlcd 

’ , ’ gammad ’ , 

* ’ rgas 

* ’dhldrd 

’,’p00d 

’ , ’mbard 

’ , ’ rbard 

’,’dcsdrd ’, 

’ , ’rholod 

’,’ulod 

’ , ’ pchmb 

’,’tchmb ’/ 


1 FORMAT (1615) 

2 FORMAT(4E15.6) 

3 FORMAT (3E 15. 6) 

4 FORMAT (A) 

5 FORMAT ( ’ Enter X (ft), P (lbf/ft~2), and T (*R) for point \ 

* 13,’ ’) 

6 FORMAT ( 1 P4E 15.6) 

7 FORMAT ( 2X , A8 , 2X , A8 , 2X , A8 , 2X , A8 , 2X , A8 ) 

8 FORMAT (2X,A8,1PE13.5,2X,A8,E13.5,2X f A8,E13.5) 

9 FORMAT ( 1P3E15. 6) 

IF(IGO.EQ.I) THEN 

WRIT£(*, ’ (A\) ’ ) ' Do you wish to use old data with or without chan 
*ges? Y or N ’ 

READ(*,4)ANS 

IF(ANS. EQ. ’Y’ .OR.ANS. EQ. ’y’ ) GO TO 24 
ENDIF 
IG0=1 

WRITE(*,*) ’ ’ 

WRITE(* , ’ (A\) ’ ) ’ Is your rocket input on file? Y OR N ’ 
READ(*,4)ANS 

IF(ANS.EQ. ’ Y’ .OR. ANS. EQ. ’ y ’ ) THEN 

WRITE(*, ’ (A\) ’ ) ' Does the file need to be rewound? Y OR N ’ 
READ(*,4)ANS 

IF(ANS. EQ. ’Y’ .OR.ANS. EQ. ’y’ ) REWIND 5 
READ(5,4,END=99)TITLE 
READ(5, 1 ,END=99)NVAL 
IF(NVAL.EQ.O) GO TO 99 

READ(5, 3)(XBARD(I) , PBAR(I) ,TBAR(I) ,1=1 ,NVAl) 

PCHMB=PBAR( 1 ) 

TCHMB=TBAR(1 ) 

READ( 5 , 2 ) ND , TAUD , DTAUD , NRD 
REA0( 5 , 2 ) LAMDAD , MUD 
READ(5,2)CDIAM,TDIAM,XLCD 
READ( 5 , 2 )GAMMAD , RGAS , POOD 
READ(5,2)MBARD,RBARD 
READ ( 5 , 2 ) DCSDRD , DHLDRD , RHOLOD , ULOD 
ELSE 

WRITE(», ’(A\)’)’ How many points along centerline? ’ 
READ(*,*,END=99)NVAL 
IF(NVAL.EQ.O) GO TO 99 
DO 21 1=1 ,NVAL 
WRITE(*,5)I 

READ ( * , * ) XBARD ( I ) , PBAR ( I ) , TBAR ( I ) 

21 CONTINUE 

PCHMB=PBAR( 1 ) 

TCHMB=TBAR(1 ) 

WRITE(*,*)’ Enter Title’ 
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READ(*,4)TITLE 

WRITE(*,*)’ Enter N (pressure interaction index) and NR’, 

* ’ (enthalpy interaction index)’ 

READ(*,*)ND,NR 

WRITE(*,*)’ Enter TAU (sensitive time lag - sec) and DTAU’ , 

* ’ (invarient time lag - sec)’ 

READ(*,*)TAUD,DTAUD 

WRITE(*,*)’ Enter LAMDA and MU (real and imaginary parts’, 

* ’ of frequency’ 

READ (*,*)LAMDAD, MUD 

WRITE(*,*)’ Enter XLCD (length of combustion chamber - ft)’ 
READ(*,*)XLCD 

WRITE(*,*)’ Enter CDIAM (chamber diameter - ft) and TDIAM’ , 

* ’ (throat diameter - ft)’ 

READ(*,*)CDIAM, TDIAM 

WRITE(*,*) ’ Enter GAMMA (ratio of specific heats), RGAS ’ , 

* ’ (gas constant - ft"2/sec~2/°R) ’ 

READ(*,*)GAMMAD,RGAS 

WRITE(*,*)’ Enter POO (maximum overpressure - lbf/ft~2)’ 

READ (*,*) POOD 

WRITE(*,*)’ Enter MBAR (mean combustion response function -’ 

* ’ Ibm/sec)’ 

WRITE(*,»)’ and RBAR (mean mixture ratio)’ 

READ ( * , * ) MBARD , RBARD 

WRITE(*,*)’ Enter DCSDR (dc*/dr - ft/sec) and DHLDR’ , 

* ’ (dh /dr - ft~2/sec^2) ’ 

READ(*,*)DCSDRD,DHLDRD 

WRITE(*,*) ’ Enter RHOLO (mass of liquid/unit chamber vol 

* ’ lbm/ft*3) ’ 

WRITE(* , *) ’ and ULO (axial component of liquid velocity’, 

* ’ - ft/sec)’ 

READ ( * , * ) RHOLOD , ULOD 
WRITE(5,4)TITLE 
WRITE(5, 1 )NVAL 

WRITE ( 5,9) (XBARD (I) , PBAR ( I ) ,TBAR( I) , 1=1 ,NVAL) 

WRITE(5, 6)ND, TAUD , DTAUD , NR 
WRITE ( 5 , 6 ) LAMDAD , MUD 
WRITE(5,6)CDIAM, TDIAM, XLCD 
WRITE(5,6)GAMMAD, RGAS, POOD 
WRITE(5,6) MBARD , RBARD 
WRITE ( 5 , 6 ) DCSDRD , DHLDRD , RHOLOD , ULOD 
ENDIF 

CALL NONDIM(HOLD) 

RETURN 
24 CONTINUE 

WRITE(*,’(A\)’)’ are there any changes? Y or N ’ 

READ(*,4)ANS 

IF(ANS. EQ. ’N’.OR.ANS.EQ. ’n’) THEN 
CALL NONDIM(HOLD) 

RETURN 

ENDIF 

WRITE(«, ’(A\)’)’ Do you wish to change title? Y or N ’ 
READ(*,4)ANS 

IF(ANS.EQ. ’Y’ .OR.ANS.EQ. ’y’) THEN 
WRITE(*,*)’ Enter Title’ 

READ(*,4)TITLE 
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ENDIF 
GO TO 29 

27 CONTINUE 

WRITE(*,*)’ VARIABLE NAMES AND DESCRIPTIONS’ 

WRITE(*,*)’ ’ 

WRITE(*,*)’ NO - pressure Interaction index’ 
WRITE(*,*)’ TAUD - sensitive time lag 
WRITE(*,*) ’ DTAUD - invariant time lag 
WRITE(*,*)’ NRD - enthalpy interaction Index’ 

WRITE(*,*) ’ LAMDAD - damping of perturbation’ 

WRITE(*,*) ’ MUD - frequency of perturbation’ 

WRITE(*,*) ’ CDIAM - chamber diameter 
WRITE(*,*)’ TDIAM - throat diameter 
WRITE(*,*)’ XLCD - length of combustion chamber 
WRITE(*,*) ’ GAMMAD - ratio of specific heats’ 

WRITE(*,*) ’ RGAS - gas constant 

* ’ ( ft/sec )" 2/ °R’ 

WRITE(*,*)’ POOD - maximum pressure 

* ’ lbf/ft~2’ 

WRITE(*,*)’ MBARD - mean combustion response funct. 

* ’ 1 bm/sec ’ 

WRITE(*,*) ’ RBARD - mean mixture ratio’ 

WRITE(*,*)’ DCSDRD - d(c*)/d (mixture ratio) 

WRITE(*,*)’ DHLDRD - d(enthalpy)/d(mixture ratio) 

* ’ft~2/sec~2’ 

WRITE(*,*) ’ RHOLOD - mass of liquid/unit chamber volume 

* ’ lbm/ft~3’ * 

WRITE(* , *) ’ ULOD - axial component of liquid velocity 
WRITE(*,*) ’ PCHMB - chamber pressure at injector 

* ’ lbf/ft~2 ’ 

WRITE(*,*)’ TCHMB - chamber temperature 
WRITEO,*)’ ' 

GO TO 30 

28 CONTINUE 

WRITE(*,*)’ VARIABLE NAMES AND VALUES’ 

WRITE(*,*) ’ ’ 

WRITE(*,8)(VAR(I),HOLD(I) ,1=1,20) 

29 CONTINUE 
WRITE(*,*) ’ ’ 

WRITE(*,*)’ Enter ? to print variable names & descriptions’ 

WRITEC*,*)’ to print variable names & values’ 

WRITE(*,*)’ END when all changes have been made’ 

WRITE(*,*) ’ ’ 

30 CONTINUE 

WRITE(*, ’ (A\) ’ ) ’ Enter variable name and new value, END, 

* ’ 

CALL ZREADCNAME, VALUE) 

IF(NAME.EQ. ’?’) GO TO 27 

IF(NAME. EQ. ’#’ ) GO TO 28 

IF(NAME.EQ. ’ END’ .OR. NAME. EQ. ’end’) THEN 
CALL NONDIM(HOLD) 

RETURN 

ENDIF 

DO 31 11=1,20 
I=II 

I F ( NAME . EQ . VARP ( I ) . OR . NAME . EQ . VARL ( I ) ) GO TO 32 


sec’ 

sec’ 


ft’ 

ft’ 

ft’ 

9 

9 

9 

9 

9 

9 


ft/sec’ 

9 

9 

9 

9 

ft/sec’ 

9 

9 

°R’ 


?, or # 
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31 CONTINUE 

WRITE(*,*)' Invalid name, try again’ 
GO TO 27 

32 CONTINUE 
HOLD(I)=VALUE 
GO TO 30 

99 CONTINUE . 

STOP 

END 
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SUBROUTINE SETVAL(VAl,ID) 

COMMON /DIMVAL/HOLDD( 20 ) , XBARD( 50 ) , PBAR( 50 ) , TBAR( 50 ) 
VAL=HOLDD(ID) 

RETURN 

END* 
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SUBROUTINE SETVAR(VAL, ID) 

COMMON /CMPVAL/X 1 , Y 1 , Z 1 , W 1 , M 1 , PO , P 1 , UO , U 1 , RFH , RFK , RFP , 

* S , GF, GOX , RFA , RFC 

COMMpN /RELVAL/N , TAU , DTAU , NR , RBAR , MBAR , GAMMA , POO , DHLDR , CSTAR , 

* DCSDR , RHOLO , ULO , LAMDA , MU , T AUT , UBAR ( 50 ) , XBAR ( 50 ) , XLC 
COMMON /RESULT/PP,UP,SIGP,FUNB 

COMMON /INTVAL/NVAL 

COMMON /DIMVAL/HOLDD ( 20 ) , XBARD ( 50 ) , PBAR ( 50 ) , TBAR ( 50 ) 

REAL MBAR , N , NR , LAMDA , MU , RVAR (13) 

REAL MBARD , ND , NRD , LAMDAD , MUD 

COMPLEX S,X1,Y1,Z1,W1,M1,P0,P1,U0,U1 ,GF, GOX, RFH, RFK, RFP, RFA, RFC 
COMPLEX FP1 , FU 1 , FSIGP , PP , UP , SIGP , FUNB , CSINH , CCOSH , CVAR (17) 
EQUIVALENCE (N,RVAR( 1 ) ) , (XI ,CVAR( 1 ) ) 

EQUIVALENCE 

* ( ND , HOLDD( 1 ) ) , (TAUD , HOLDD( 2 ) ) , ( DTAUD , HOLDD ( 3 ) ) , 

* (NRD,H0LDD(4)) , (LAMDAD, HOLDD(5)) , (MUD,HOLDD(6)) , 

* (CDIAM,HOLDD(7)) , (TDIAM,HOLDD(8) ) , (XLCD,HOLDD(9)) , 

* (GAMMAD , HOLDD( 1 0 ) ) , ( RGAS , HOLDD( 1 1 ) ) , ( POOD , HOLDD ( 12)), 

* ( MBARD , HOL DD ( 1 3 ) ) , ( RBARD , HOL DD ( 1 4 ) ) , ( DCSDRD , HOL DD ( 1 5 ) ) , 

* (DHLDRD,HOLDD( 16)), (RHOLOD,HOLDD( 17)) , (ULOD,HOLDD(18)) , 

* (PCHMB,HOLDD(19) ) , (TCHMB,HOLDD(20)) 

DATA PI/3. 141593/, GC/32. 174/ 

HOLDD(ID)=VAL 

IF(ID.EQ.I) THEN 
C ND 

N=ND 
RETURN 
ENDIF 

IF(ID.EQ.2) THEN 
C TAUD 

ASTAR=SQRT (GAMMAD*RGAS*TBAR( 1 ) ) 

TAU=TAUD*ASTAR/XLCD 

TAUT=TAU+DTAU 

RETURN 

ENDIF 

IF(ID.EQ.3) THEN 
C DTAUD 

ASTAR=SQRT(GAMMAD*RGAS*TBAR(1 )) 

DT AU= DT AUD*AST AR/XLCD 
TAUT=TAU+DTAU 
RETURN 
ENDIF 

IFUD.EQ.4) THEN 
C NRD 

NR=NRD 
RETURN 
ENDIF 

IF(ID.EQ.5) THEN 
C LAMDAD 

ASTAR=SQRT(GAMMAD*RGAS*TBAR(1 )) 

LAMDA= LAMDAD*XLCD/ASTAR 
S=CMPLX( LAMDA, MU) 

RETURN 

ENDIF 

IFUD.EQ.6) THEN 
C MUD 
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ASTAR=SQRT(GAMMAD*RGAS*TBAR( 1 ) ) 

MU=MUD*XLCD*PI/ASTAR 

S= CMPI.X ( LAMDA , MU ) 

RETURN 

ENDfF 

IF( ID. EQ. 7) THEN 
C CDIAM 

CAREA=0. 25*PI*CDIAM**2 
ASTAR=SQRT(GAMMAD*RGAS*TBAR( 1 ) ) 

DO 21 1=1 ,NVAL 

RHOBAR=PBAR ( I ) *GC/ ( RGAS*T BAR ( I ) ) 
UBARD=MBARD/ ( RHOBAR*CAREA ) 

UBAR ( I ) =UBARD/AST AR 

21 CONTINUE 

RETURN 

ENDIF 

IFUD.EQ.8) THEN 
C TDIAH 

TAR£A=0. 25*PI*TDIAM**2 
ASTAR=SQRT(GAMMAD*RGAS*TBAR( 1 )) 
CSTARD= PBAR (NVAL) *T AREA*GC/MBARD 
CST AR= CST ARD/AST AR 
RETURN 
ENDIF 

IF(ID.EQ.9) THEN 
C XLCD 

ASTAR=SQRT(GAMMAD*RGAS*TBAR( 1 ) ) 

TAU=TAUD*ASTAR/XLCD 

DTAU= DTAUD*ASTAR/XL CD 

TAUT=TAU+DTAU 

LAMDA=LAMDAD*XLCD/ASTAR 

MU=MUD*XLCD*PI/AST AR 

S=CMPLX( LAMDA, MU) 

DO 22 1=1, NVAL 
XBAR(I)=XBARD(I )/XLCD 

22 CONTINUE 

RETURN 

ENDIF 

IF(ID.EQ.IO) THEN 
C GAMMAD 

GAMMA= GAMMA D 
CAREA=0. 25*PI*CDIAM**2 
TAREA=0. 25*PI*TDIAM*<2 
ASTAR=SQRT(GAMMAD*RGAS*TBAR( 1 )) 

T AU=TAUD*AST AR/XLCD 

DTAU= DTAUD* AST AR/XLCD 

TAUT=TAU+DTAU 

LAMDA= LAMDAD*XLCD/AST AR 

MU=MUD*XLCD*PI/ASTAR 

S=CMPLX( LAMDA, MU) 

ULO=ULOD/ASTAR 

DCSDR= DCSDRD/ AST AR 

RHOB 1 = PBAR ( 1 ) *GC/ ( RGAS*T BAR ( 1 ) ) 

MBAR=MBARD/(RHOB1*ASTAR*CAREA/XLCD) 

CSTARD= PBAR ( NVAL ) *T AREA*GC/MBARD 

CST AR= CST ARD/AST AR 
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DO 23 1=1 , NVAL 

RHOBAR=PBAR( I)*GC/(RGAS*TBAR( I ) ) 
UBARD=MBARD/ ( RHOBAR*CAREA ) 
UBAR(I)=UBARD/ASTAR 

23 CONTINUE 

RETURN 

ENDIF 

IF(ID.EQ.II) THEN 
C RGAS 

CAREA=0. 25*PI*CDIAM**2 
TAREA=0. 25*PI*TDIAM**2 
ASTAR=SQRT(GAMMAD*RGAS*TBAR( 1 ) ) 

T AU=TAUD*AST AR/XLCO 
DT AU= DT AUD*AST AR/X LCD 
TAUT=TAU+DTAU 
L AMDA= LAMDAD*XLCD/ AST AR 
MU=MUD*XLCD*PI/ASTAR 
S=CMPLX( LAMDA , MU) 

ULO=ULOD/ASTAR 
OCSDR= OCSDRD/ AST AR 
RHOB1=PBAR(1 )*GC/(RGAS*TBAR(1 )) 
RHOLO= RHOLOD/RHOB 1 
MBAR=MBARD/(RHOB1*ASTAR*CAREA/XLCD) 
CSTARD= PBAR ( NVAL ) *TAREA*GC/MBARD 
CST AR=CST ARD/ AST AR 
DO 24 1=1, NVAL 

RHOBAR=PBAR( I )*GC/(RGAS*TBAR( I ) ) 
UBARD=MBARD/(RHOBAR*CAREA) 

UBAR ( I ) = UBARD/ AST AR 

24 CONTINUE 

RETURN 

ENDIF 

IF(ID.EQ.12) THEN 
C POOD 

POO=POOD/PCHMB 
RETURN 
ENDIF 

IF(ID. EQ. 13) THEN 
C MBARD 

CAREA=0. 25*PI*CDIAM**2 
TAREA=0. 25*PI*TDIAM**2 
ASTAR=SQRT(GAMMAD*RGAS*TBAR( 1 ) ) 

RHOB 1 = PBAR ( 1 ) *GC/ ( RGAS*TBAR ( 1 ) ) 
MBAR=MBARD/ ( RHOB1 *ASTAR*CAREA/XLCD ) 
CSTARD= PBAR ( NVAL ) »TAREA*GC/MBARD 
CST AR= CST ARD/AST AR 
DO 25 1=1, NVAL 

RHOBAR= PBAR ( I ) *GC/ ( RGAS*TBAR (I)) 
UBARD=MBARD/ ( RHOBAR*CAREA ) 

UBAR ( I )= UBAR D/AST AR 

25 CONTINUE 

RETURN 

ENDIF 

IF(ID. EQ. 14) THEN 
C RBARD 

RBAR=RBARD 
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RETURN 

ENDIF 

IF(ID. EQ. 15) THEN 
C DCSDRD 

AST>R=SQRT(GAMMAD*RGAS*TBAR( 1 ) ) 
DCSDR= DCSDRD/ AST AR 
RETURN 
ENDIF 

IF(ID. EQ. 16) THEN 
C DHLDRD 

DHLDR=DHLDRD 
RETURN 
ENDIF 

IF(ID. EQ. 17) THEN 
C RHOLOD 

RHOB 1 = PBAR ( 1 ) *GC/ ( RGAS*TBAR ( 1 ) ) 
RHOLO=RHOLOD/RHOB1 
RETURN 
ENDIF 

IFUD.EQ.18) THEN 
C ULOD 

ASTAR=SQRT(GAMMAD*RGAS*TBAR( 1 ) ) 
ULO=ULOD/ASTAR 
RETURN 
ENDIF 

IF(ID.EQ.19) THEN 
C PCHMB 

CAREA=0. 25*PI*CDIAM**2 
TAREA=0. 25*PI*TDIAM«*2 
ASTAR=SQRT(GAMMAD*RGAS*TBAR( 1 ) ) 
FAC=PCHMB/PBAR( 1 ) 

DO 26 1= 1 , NVAL 
PBAR(I)=FAC*PBAR(I) 

RHOBAR= PBAR ( I ) *GC/ ( RGAS*T BAR (I)) 
UBARD=MBARD/(RHOBAR*CAREA ) 

UBAR ( I ) =UBARD/AST AR 

26 CONTINUE 

CSTARD= PBAR ( NVAL ) *TAREA*GC/MBARD 
CST AR= CST ARD/AST AR 
RHOB 1 = PBAR ( 1 ) *GC/ ( RGAS*TBAR ( 1 ) ) 
RHOLO= RHOLOD/RHOB 1 
MBAR=MBARD/ ( RHOB 1 *ASTAR*CAREA/XLCD) 
POO=POOD/PCHMB 
RETURN 
ENDIF 

IF(ID. EQ. 20) THEN 
C TCHMB 

DO 27 1=1, NVAL 
TBAR ( I ) = FAC*TBAR ( I ) 

27 CONTINUE 

CAREA=0. 25*PI*CDIAM**2 
TAREA=0. 25*PI*TDIAM**2 
ASTAR=SQRT(GAMHAD*RGAS*TBAR( 1 ) ) 
FAC=TCHMB/TBAR( 1 ) 

DO 28 1=1, NVAL 

RHOBAR= PBAR (I)*GC/( RGAS*TB AR (I)) 
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UBARD=MBARD/ ( RHOBAR*CAREA ) 

UBAR ( I ) =UBARD/AST AR 
28 CONTINUE 

CSTARD= PBAR ( NVAL ) *TAREA*GC/MBARD 
CSTAR= CST ARD/AST AR 
RHOB 1 = PBAR ( 1 ) *GC/ ( RGAS*T BAR ( 1 ) ) 
RHOLO= RHOLOD/RHOB 1 
MBAR=MBARD/(RHOB1*ASTAR*CAREA/XLCD) 
ENDIF 
RETURN 
END 


/ 
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SUBROUTINE ZREAD(NAME, VALUE) 

SUBROUTINE TO READ NAMES AND VALUES 
CHARACTER* 1 NAME(8) 

CHARACTER* 1 CARD(80) , PLUS, MINUS , PERIOD, LE, E, NUMBER ( 10) 

CHARACT ER*1 LEND(3),CEND(3), POUND , QUEST , BL K , COMMA 
CHARACTER*80 DCARD 
EQUIVALENCE ( CARD (1), DCARD) 

DATA PLUS/ ’ + V , MINUS/ ’ - ’ / , PERIOD/ ’ . ’/ , LE/ ’ e ’ / , E/ ’ E ’ / , BLK/ ’ 7 
DATA NUMBER/ , 0 , t , r, , 2 , , , 3 , f , 4 , f , 5 , ,’6* t , 7 , , , 8 , , , 9 , /,C0MMA/ , l V 
DATA LEND/’e’, ’n’, ’d’/.CEND/’E’ , ’N’.’D’/, POUND/’#’/, QUEST/’?’/ 

1 FORMAT (A) 

DO 21 1=1,8 
NAME(I)=BLK 
21 CONTINUE 

READ C* , 1 ) DCARD 
IF(CARD(1).EQ. POUND) THEN 
NAME (IMPOUND 
RETURN 


ENDIF 

IF(CARD(1).EQ. QUEST) THEN 
NAME(1)=QUEST 
RETURN 
ENDIF 

DO 22 1=1,3 

IF(CARD(I).NE.LEND(I).AND.CARD(I).NE.CEND(I)) GO TO 23 
NAME(I)=CEND(I) 

22 CONTINUE 
RETURN 

23 CONTINUE 

DO 24 1=1,8 
II=I 

IF(CARD(I).EQ.BLK. OR. CARD(I).EQ. COMMA) GO TO 25 
NAME(I)=CARD(I) 

24 CONTINUE 

25 CONTINUE 

DO 26 1=11,80 
ID=I 

IF(CARD(I).NE.BLK. AND. CARD(I).NE. COMMA) GO TO 27 

26 CONTINUE 
VALUE=0.0 

WRITE(*,*)’ NO VALUE GIVEN, ZERO ASSUMED’ 

RETURN 

27 CONTINUE 
SIGN=1 . 0 

IF(CARD(ID).EQ. MINUS) THEN 
SIGN=-1 .0 


ID=ID+1 

ELSEIF(CARD(ID) . EQ. PLUS) THEN 
ID=ID+1 
ENDIF 


WHOLE=0.0 
DO 30 I=ID,80 
II=I 

IF(CARD(I).EQ. PERIOD) GO TO 31 
IF(CARD(I ) . EQ. PLUS) GO TO 36 
IF(CARD(I).EQ. MINUS) GO TO 36 
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I F ( CARO ( I ) . EQ. E.OR.CARD(I) -EQ. LE) GO TO 35 
DO 28 J=1 , 10 
JJ=J-1 

IF(CARD(I).EQ.NUMBER(J)) GO TO 29 

28 CONTINUE 

VALU£=SIGN*WHOLE 

IF(CARD(I).EQ.BLK) RETURN 

WRITE(*,*)’ INPUT ERROR, VALUE SET TO ZERO’ 

VALUE=0.0 

RETURN 

29 CONTINUE 

WHOLE=WHOLE*10. 0+JJ 

30 CONTINUE 
VALUE=SIGN*WHOLE 
RETURN 

31 CONTINUE 
ID=II+1 
FRACT=0.0 
I COUNT =0 

DO 34 I=ID,80 
I COUNT = I COUNT + 1 
II=I 

IF(CARD(I).EQ. PERIOD) THEN 
WRITE(*,*) ’ INPUT ERROR, VALUE SET TO ZERO’ 
VALUE=0.0 
RETURN 
ENDIF 

IF(CARD(I) . EQ.PLUS) GO TO 36 
IF(CARD(I).EQ. MINUS) GO TO 36 
IF(CAR0(I) . EQ, E.OR. CARD(I) . EQ. LE) GO TO 35 
DO 32 J=1,10 
JJ=J-1 

IF(CARD(I).EQ.NUMBER(J)) GO TO 33 

32 CONTINUE 

VALUE=SIGN* (WHOLE+FRACT) 

IF(CARD(I).EQ.BLK) RETURN 

WRITE(*,*)’ INPUT ERROR, VALUE SET TO ZERO’ 

VALUE=0. 0 

RETURN 

33 CONTINUE 

FRACT=FRACT+JJ/10.0**ICOUNT 

34 CONTINUE 

VALUE=SIGN* (WHOLE+FRACT ) 

RETURN 

35 CONTINUE 
11=11+1 

36 CONTINUE 

VALUE=SIGN* (WHOLE+FRACT) 

SIGN=1.0 

IF(CARD(II).EQ. MINUS) THEN 
SIGN=-1 .0 
11=11+1 

ELSEIF(CARD(II). EQ.PLUS) THEN 
11=11+1 
ENDIF 
WHOLE=0.0 
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DO 39 1=11,80 
DO 37 J=1,10 
JJ=J-1 

IF(CARD(I).EQ.NUMBER(J)) GO TO 38 

37 CONTINUE 

VALUE=VALUE*10.0**(SIGN*WHOLE) 

IF(CARD(I) . EQ.BLK) RETURN 

WRITE(*,*)' INPUT ERROR, VALUE SET TO ZERO' 

VALUE=0.0 

RETURN 

38 CONTINUE 

WHOLE=WHOLE* 1 0 . 0+ J J 

39 CONTINUE 

VALUE=VALUE* 1 0 . 0** (SIGN*WHOLE) 

RETURN 

END 
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APPENDIX V 

FORTRAN Listing of FEEDLINE 
Wilbur C. Armstrong 
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Listing of Program FEEDLINE 


PROGRAM TO CALCULATE ADMITTANCE COEFFICIENTS FROM TANK 
THROUGH FEEDLINES TO ROCKET CHAMBER 


A = VELOCITY OF SOUND IN THE FLUID (FT/S) 

AREA(I) = AREA OF FEED LINE (FT"2) 

CMAN = CAPACITANCE ASSOCIATED WITH THE MANIFOLD 
CTANK = CAPACITANCE ASSOCIATED WITH THE FUEL TANK 
DENS = DENSITY OF FLUID (LB/FT~3) 

DIA(I) = DIAMETER OF LINE 
DPROR = PRESSURE DROP ACROSS ORIFICES (LB/FT~2) 

G ( I ) = OVERALL ADMITTANCE AS SEEN FROM ENGINE 
GRAV = ACCELERATION OF GRAVITY (32.2 FT/S~2) 

INERT = INERTANCE DUE TO DUCT CURVATURE 
KMAN = BULK MODULUS OF FLUID AT MANIFOLD PRESSURE (LB/FT~2) 
KTANK = BULK MODULUS OF FLUID IN TANK (LB/FT~2) 

L ( I ) = LENGTH OF LINE (FT) 

LFLOW = FLOW RATE OF FLUID THROUGH THE LINE (LB/S) 

PCHMB = PRESSURE INSIDE COMBUSTION CHAMBER (LB/FT2) 

RATIO = RATIO OF INNER RADIUS TO OUTER RADIUS 
TFLOW = TOTAL FLOW RATE OF FLUID INSIDE ENGINE(LB/S) 

TL = TIME CONSTANT 
VOL = VOLUME OF FUEL TANK (FT~3) 

VOLMF = VOLUME OF MANIFOLD (FT^3) 

ZLINE = IMPEDANCE ASSOCIATED WITH THE FEED LINE 
ZOR = IMPEDANCE ASSOCIATED WITH THE ORIFICES 


COMPLEX G ( 76 ) ,CTANH,G1 ,S 

REAL AREA (75) ,DIA(75) , L(75) , INERT, INRAD, KMAN, KTANK, LBEND, LFLOW, 

* LFREQ,LPRME, MAG, NEWLN, MAGI 

INTEGER N ( 75 ) ,SECTN(75) ,PTS,RSPON,SECT,SEGMN 
CHARACTER ANS*1 

CHARACTER*20 TITLE, TITLF,TITLO 
DATA GRAV/32. 2/, PI/3. 1415927/ 

DATA TITLF/’ FUEL LINE 7 

DATA TITLO/’ OXYGEN LINE 7 

1 FORMAT ( E 1 5 . 6 ) 

2 FORMAT ( 15 , 3E 1 5.6) 

3 FORMAT ( 1 P4E 1 5 . 6 ) 

4 FORMAT ( 1 PE 1 3 . 5 , ’ ( ' , E12 . 5, ' , ’ , E12. 5, ' ) ( ’ , E12.5, ' , ’ , E12. 5 , ’ ) ’ ) 

5 FORMAT (/’ FREQ’ ,8X, ’FREQ-NORM’ , 9X,’G(R)’,11X,’G(I)’/) 

6 FORMAT ( /2x, FREQ" ’ ,7X, ’"FREQ-NORM" ’, 5X, ’ " /G1/"’,6X, 

* ’" /G/"’/) 

7 FORMAT ( ’ " ’ , A , ”” ) 

8 FORMAT ( 1 PE 15.6) 

9 FORMAT (15, 1P3E15.6) 

WRITE(*,’(A\)’)’ IS THIS SETUP FOR FUEL OR OXIDIZER? ENTER F OR 0 

*. ’ 

READ(», ’ (A) ’ )ANS 

IF(ANS.EQ. ’F’.OR.ANS.EQ.’f’) THEN 
TITLE=TITLF 

OPEN(UNIT=1 1 , FILE= ’ FUEL.DAT ’ ,STATUS= ’OLD’ ) 

ELSE 

TITLE=TITLO 
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fc 

L 


0PEN(UNIT=1 1 , FILE= ’ LOX. DAT ’ , STATUS: ’OLD’ ) 

ENDIF ' 

OPEft( UNIT= 1 2 , FI LE= ’ ADMIT . DAT ’ , STATUS= ’ NEW ’ ) 

OPEN ( UNIT= 7 , FILE= * ADMIT . PRN ’ , STATUS= ’ NEW ’ ) 

WRITE(*, ’ (A\) ’ ) ’ IF THERE IS DATA STORED ENTER YES ’ 

READ(*, ’ (A) ’ )ANS 

IF(ANS.EQ. ’N’.OR.ANS.EQ. ’n’) THEN 
RSPON=3 
GO TO 21 
ENDIF 

READ(11, 1)DENS 
READ( 11,1 )TFLOW 
READ(1 1 , 1 )VOLMF 
READ( 11,1 )KMAN 
READ( 11,1 )PCHMB 
READ( 11,1 )DPROR 
READ(11 ,1)VOL 
READ(11, 1)LFLOW 
READd 1 , 1 )KTANK 
READ(11, 1)A 
READ(1 1,1)CTANK 
READ(1 1,1) OMAN 
READ(1 1 ,2)SEGMN 

READ( 1 1 ,2) (SECTN(I ) , L(I) ,AREA(I) ,DIA(I) ,1=1 , SEGMN) 

The first stage in this program is to define the parameters then 
we will begin the initial calculations. Because these parameters 
are as likely to change as not, a provision is made to update the 
parameters if necessary. 

WRITE(12,*) ’ * 

WRITE(12,*)TITLE 
WRITE(12,*)’ ’ 

WRITE(1 2,*)’ PRESENT CONDITIONS ARE AS FOLLOWS:’ 

WRITE (1 2 ,*)’ DENS= DENS 

WRITEd2,*) ’FUEL TANK VOLUME=’,VOL 

WRITE (1 2 , * ) ’ MANIFOLD VOLUME= * , VOLMF 

WRITE (12,*) ’ENGINE CHAMBER PRESSURE: ’ , PCHMB 

WRITE (12,*) ’PRESSURE DROP ACROSS ORIFICE= ’ ,DPROR 

WRITE (12,*) ’TOTAL FLOW RATE=’,TFLOW 

WRITE(12,*) ’LINE FLOW RATE=’,LFLOW 

WRITE (12,*) ’BULK MOD. OF FUEL TANK=’ , KTANK 

WRITE(12,*) ’BULK MOD. OF MANIFOLD: ’, KMAN 

WRITE( 12,*) ’VELOCITY OF SOUND IN FLUID:’, A 

WRITE (12,*) ’CAPACITANCE OF FUEL TANK:’ ,CTANK 

WRITE(12,*) ’CAPACITANCE OF MANIFOLD: ’ ,CMAN 

WRITE (12,*)’ STATUS LENGTH AREA DIAMETER’ 

WRITE(12 , 9)(SECTN(I) ,L(I) ,AREA(I) ,DIA(I) ,1=1 ,SEGMN) 

WRITE(12,*)’ ’ 

WRITE(*,*) ’ ’ 

WRITE(*,*)TITLE 
WRITE(*,*) ’ ’ 

WRITE(*,») ’PRESENT CONDITIONS ARE AS FOLLOWS:’ 

WRITE(*,*) ’DENS=’, DENS 
WRITE(*,*) ’FUEL TANK VOLUME:’, VOL 
WRITE (*,*) ’MANIFOLD VOLUME: ’, VOLMF 
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WRITEC*,*) ’ENGINE CHAMBER PRESSURE= ’ ,PCHMB 

WRITEC*,*) ’PRESSURE DROP ACROSS ORIFICE= ’ , DPROR 

WRITE (*,*) ’TOTAL FLOW RATE=’,TFLOW 

WRITEC*,*) ’LINE FLOW RATE=’,LFLOW 

WRITEC* , *) ’BULK MOD. OF FUEL TANK= ’ , KTANK 

WRITEC*,*) ’BULK MOD. OF MANIFOLD= ’ , KMAN 

WRITEC*,*) ’VELOCITY OF SOUND IN FLUID=’,A 

WRITEC*,*) ’CAPACITANCE OF FUEL TANK=’,CTANK 

WRITEC*,*) ’CAPACITANCE OF MANIFOLD= ’ ,CMAN 

WRITEC*,*)’ STATUS LENGTH AREA DIAMETER’ 

WRITEC*, 9) (SECTN(I), L(I),AREA(I),DIA(I), 1=1, SEGMN) 

WRITEC*,*)’ If revisions on the design have been made’ 

WRITEC*,*) ’ (changes in fuel, pipe length, diameter, bends, etc.)’ 
WRITEC*, ’ (A\) ’ ) ’ Please enter yes for revisions or no to continue. 
* » 

READC*, ’ (A) ’ )ANS 

IF(ANS.EQ. ’N’.OR.ANS.EQ. ’n’) GO TO 28 

WRITE (*,*) ’CHANGES IN FUEL LINE SPECIFICATIONS ENTER 1’ 

WRITEC*,*) ’CHANGES IN ENGINE, TANK, AND FUEL SPECS. ENTER 2’ 
WRITEC*,*) ’CHANGES FOR BOTH FUEL LINE AND ENGINE SPECS. ENTER 3’ 
READ(*,*)RSPON 
REWIND 11 

IFCRSPON.EQ.2) GO TO 27 
C 

C AT THIS POINT, YOU NEED TO FURNISH ME WITH SOME INFORMATION, LIKE 
C WHAT ARE THE LENGTHS AND AREAS OF THE SEGMENTS, HOW MANY BENDS ARE 

C THERE, WHERE ARE THEY LOCATED ON THE PIPE, WHAT IS THE ANGLE THAT 

C THEY MAINTAIN, WHAT IS THE RADIUS FOR THE BEND, AND ARE THERE 
C VALVES, WHERE ARE THEY LOCATED, ETC. 

C 

21 CONTINUE 

*** REMEMBER TO RESOLVE AMBIGUITY CONCERNING THE LENGTHS FOR 
MANIFOLDS AND ORIFICES, SINCE THE PARAMETERS FOR THEIR 
IMPEDENCE WILL BE INPUT. 2-22-90 

WRITE(* , ’ (A\) ’ ) ’HOW MANY SEGMENTS IS THE PIPE BROKEN INTO? ’ 
READ(*,*)SEGMN 
L(I)=0.0 
AREA(I)=0.0 
DIA(I)=0.0 
SECTN(I)=0 
DO 26 1=1 ,SEGMN 

WRITEC*,*) ’SPECIFY 0 for bend, 1 FOR STRAIGHT PIPE,' 

WRITEC*,*)’ 2 FOR MANIFOLD, OR 3 FOR ORIFICE’ 

WRITE(*,*)’IN ORDER, SPECIFY LENGTH, AREA, AND DIAMETER OF’ 
WRITEC*,*) ’EACH SEGMENT’ 

WRITEC*,*) ’EX. INPUT: 1 .LENGTH, AREA, DIAMETER’ 

WRITEC*,*)’ OUTPUT: 1 LENGTH AREA DIAMETER' 

READC*,*) SECT, VALUE, AREAB, DIME 
IF(SECT.EQ.O) GO TO 24 
IF(SECT.EQ.I) GO TO 22 
IFCSECT.EQ.2) THEN 
WRITEC*, ’ (A\) ’ ) ’DENSITY OF FLUID ’ 

READC*,*) DENS 

WRITEC*, ’ (A\) ’ ) ’TOTAL FLOW RATE ’ 
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READ(*,*)TFLOW 

WRITE(*, ’ (A\) ’ ) ’VOLUME FLUID IN MANIFOLD ’ 

R5AD(*,*)VOLMF 

WRITE(* , ’ (A\) ’ ) ’ BULK MODULUS OF FLUID IN MANIFOLD ’ 

READ(*,»)KMAN 

WRITE(*, ’ (A\) ’ ) ’ENGINE CHAMBER PRESSURE ’ 

READ(*,*)PCHMB 

CMAN= ( DENS*VOLMF*PCHMB ) / ( KMAN*T FLOW ) 

WRITE(* , * ) ’ MANIFOLD CAPACITANCE= ’ , CMAN 

VALUE=0.0 

AREAB=0.0 

DIME=0.0 

SECT=2 

GO TO 23 

ELSEIF(SECT. EQ. 3) THEN 

WRITE(* , ’ (A\) ’ ) ’ PRESSURE DROP ACROSS ORIFICE ’ 

READ(*,*)DPROR 
GO TO 22 
ENDIF 

22 CONTINUE 

L(I)=VALUE 
AREA(I)=AREAB 
SECTN(I)=SECT 
DIA( I )=DIME 

23 CONTINUE 

WRITE(*,*)SECTN(I),L(I),AREA(I),DIA(I) 

GO TO 26 

24 CONTINUE 

THIS PART OF THE PROGRAM FINDS NEW LENGTHS AND AREAS FOR 
SECTIONS OF THE PIPE WITH BENDS IN THEM BY MAKING USE OF 
THE IDEA THAT A BEND CAN BE MODELED AS A STRAIGHT PIPE WITH 
DECREASED AREA AND INCREASED LENGTH. THIS SECTION LOOKS AT 
TWO CASES: MITERED AND ELBOW BENDS. 

WRITEC*, ’(A\)’)’ PLEASE ENTER LENGTH OF BEND (E.G. 1/3 LENGTH) ’ 
READ(*,*)LBEND 

WRITE (*,*) ’THE INNER AND OUTER RADIUS OF THE BEND FROM THE LOCUS’ 
WRITE(*,») ’ AN INNER RADIUS OF 0.0 DENOTES A MITERED BEND’ 

READ ( * , * ) INRAD , OTRAD 

WRITE(*,’(A\)’) ’ENTER THE INITIAL AREA OF THE PIPE ’ 
READ(*,*)ARBND 

WRITE (*,*) LBEND , INRAD , ’ FT ’ , OTRAD , ’ FT ’ , ARBND , ’ FT~ 2 ’ 

WRITE(* , ’ (A\) ’ ) ’ ENTER THE ANGLE OF THE BEND ’ 

READ(*,*)BEND 
RAT 10= INRAD/ OTRAD 
X=RATIO 

WRITE(*,*)BEND, ’DEGREES’ 

CALL GINERT(BEND,X,Y) 

25 CONTINUE 

INERT=(Y*(OTRAD-INRAD))/ARBND 

LPRME=LBEND/ARBND 

NEWLN=LPRME+INERT 

GAMMA=NEWLN/LPRME 

L ( I ) = GAMMA* LBEND 

AREA ( I ) = ARBND/SQRT ( GAMMA ) 
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26 CONTINUE 

27 CONTINUE 

IF(RSPON. EQ. 1 ) GO TO 28 
IF(*RSP0N.EQ.2) THEN 
WRITE ( * , ’ (A\) ’ ) ’ ENTER FUEL DENSITY ’ 

READ(*,*)DENS 

WRITE(*, ’ (A\) ’ ) ’ENTER TOTAL FLOW RATE INSIDE ENGINE ’ 

READ(*,*)TFLOW 

WRITE(», ’ (A\) ’ ) ’ENTER MANIFOLD VOLUME ’ 

READ(*,*)VOLMF 

WRITEC*, ’ (A\) ’ ) ’ENTER BULK MODULUS OF FLUID INSIDE MANIFOLD ’ 
READ(*,*)KMAN 

WRITEC*, ’(A\) ’)’ ENTER CHAMBER PRESSURE IN ENGINE ’ 

READ(*,*)PCHMB 

WRITEC*, ’ (A\) * ) ’ ENTER PRESSURE DROP ACROSS ORIFICE ’ 

READ(*,*)DPROR 

ENDIF 

WRITEC*, ’ (A\) ’ ) ’ENTER FUEL TANK VOLUME ’ 

READ(*,*)VOL 

WRITEC*, ’ (A\) ’ ) ’ ENTER FLOW RATE INSIDE LINE ’ 

READ(*,*)LFLOW 

WRITEC*, ’ (A\)’ VENTER BULK MODULUS OF FLUID INSIDE TANK ’ 
READC*,*)KTANK 

WRITEC*, ’ (A\) ’ ) ’ENTER VELOCITY OF SOUND IN FLUID ’ 

READ(*,*)A 

CTANK=(DENS»VOL*PCHMB)/(KTANK*TFLOW) 

CMAN=(DENS*VOLMF*PCHMB)/(KMAN*TFLOW) 

WRITE( 1 1 , 8)DENS 
WRITE( 1 1 ,8)TFLOW 
WRITE( 1 1 , 8)VOLMF 
WRITE( 1 1 , 8)KMAN 
WRITE(11,8) PCHMB 
WRITE(11 ,8)DPROR 
WRITE( 1 1 , 8)VOL 
WRITE(11 ,8)LFL0W 
WRITE(1 1 ,8)KTANK 
WRITE( 1 1 ,8)A 
WRITE(11 ,8)CTANK 
WRITE(11 ,8)CMAN 
WRITE(11,9)SEGMN 

WRITE(1 1 ,9)SECTN(I) , L(I) ,AREA(I) , DIA(I) 

THIS SECTION COMPUTES THE NEW ADMITTANCE OVER VARYING FREQUENCIES. 

28 CONTINUE 

WRITE(*,») ’ENTER RANGE OF FREQUENCIES EX: 1,2,10’ 

WRITE(*,*) ’LOW FREQ= 1 HIGH FREQ=2 #PTS=10’ 

READ(*,*)LFREQ,HFREQ,PTS 
IF(PTS.LE.I) GO TO 32 

THIS SECTION WILL COMPUTE THE ADMITTANCE RATIO FOR THE FUEL TANK 
AND THEN IT WILL COMPUTE THE ADMITTANCE RATIOS FOR EACH SEGMENT, 
SINCE THERE ARE L(I) I=1,SEGMN LENGTHS, THEN THERE WILL BE AT LEAST 
AS MANY ADMITTANCE RATIOS, THEREFORE I AM SETTING UP AN ARRAY FOR 
EACH LENGTH L(I) HAVING AN ADMITTANCE RATIO G(I). 
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SSIZE=(HFREQ-LFREQ)/(PTS-1) 

WRITE(*,*) ’STEP SIZE’ ,SSIZE 
WRITE(*,*) * * 

WRITE(*,*) ’ FREQUENCY G G1 ’ 

WRIft(*,*) * ’ 

WRITE(12,5) 

WRITE(7 ,7)TITLE 
WRITE(7 , 6) 

ZTOP=A*TFLOW/(GRAV*PCHMB) 

Z0R=2 . 0*DPR0R*TFL0W/(LFL0W*PCHMB) 

TLT=0.0 

DO 29 1=1 ,SEGMN 
TLT=TLT+L(I) 

29 CONTINUE 
TLT=TLT/(PI*A) 

DO 31 K=1 , PTS 

W=LFREQ+SSIZE*(K-1) 

S=CMPLX(0.0,W) 

G ( 1 )=CTANK*S 
G1=G(1 )+1 .0 
DO 30 1=2 ,SEGMN+1 
IF(SECTN(I-1 ) - EQ, 1 ) THEN 
ZLINE=ZT0P/AREA(I-1) 

TL= L ( I— 1 )/A 

G(I)=(1.0+CTANH(S*TL)/(G(I-1)*ZLINE))/(1.0+G(I-1)*ZLINE* 

* CTANH(S*TL)) 

ELSEIF(SECTN(I-1 ) . EQ. 2) THEN 
G(I )= 1 . 0+CMAN*S/G(I-1 ) 

ELSEIF(SECTN(I-1 ) . EQ. 3) THEN 
G ( I ) = 1 . 0/ ( 1 . 0+ZOR*G(I-1 )) 

ENDIF 

G1=G1*G(I) 

G(I)=G(I)*G(I-1 ) 

30 CONTINUE 

MAG=CABS(G(SEGMN+1 ) ) 

MAG1=CABS(G1 ) 

WN=W*TLT 

WRITE(12,3)W,WN,G(SEGMN+1 ) 

WRITE(7,3)W,WN, MAGI, MAG 
WRITE(*,4)W,G(SEGMN+1 ) ,G1 

31 CONTINUE 

32 CONTINUE 

WRITE(*,’(A\)’)’TO CONTINUE ENTER YES ’ 

READ(*, ’(A)’)ANS 

IF(ANS. EQ. ’ Y’ .OR. ANS. EQ. ’y ’ ) GO TO 28 

STOP 

END 
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COMPLEX FUNCTION CTANH(S) 
COMPLEX CCOSH,CSINH,S 
CT ANH= CSINH ( S ) /CCOSH ( S ) 
RETURN • 

END - 


COMPLEX FUNCTION CSINH(S) 
COMPLEX S 
REAL LAMDA, MU 
LAMDA=REAL(S) 

MU=AIMAG(S) 

SINHR=SINH( LAMDA )*COS (MU) 

SINHI=COSH(LAMDA)*SIN(MU) 

CSINH=CMPLX(SINHR , SINHI ) 

RETURN 

END 


COMPLEX FUNCTION CCOSH(S) 
COMPLEX S 
REAL LAMDA, MU 
LAMDA=REAL(S) 

MU=AIMAG(S) 

COSHR= COSH ( LAMDA ) *COS ( MU ) 

COSHI=SINH(LAMDA)*SIN(MU) 

CCO$H=CMPLX ( COSHR , COSHI ) 

RETURN 

END 


SUBROUTINE GINERT(BEND,X, Y) 

AO,A1,A2 = INERTANCE FIT FOR GIVEN BEND ANGLE 
B(I) = COEFFICIENT ARRAY OF INERTANCE FIT 
BEND = ANGLE OF BEND (DEGREES) 

X = RATIO OF INNER RADIUS TO OUTER RADIUS (RI/RO) 
Y = INERTANCE 

DIMENSION B(3) 

DATA B/0 . 0 , 0 . 78770 1 4E-02 ,-0.281 4679E-04/ 

A=B( 1 )+(B(2)+B(3)*BEND)*BEND 

Y=A*(X-1.0)**2 

RETURN 

END 
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